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1.  SUMMARY 


Development  of  high-resolution  seismic  velocity  and  attenuation  models  in  the 
Qinghai-Tibetan  plateau  and  adjacent  regions  of  western  China  are  critical  to  monitoring 
efforts  because  of  their  close  proximity  (~380  km)  to  the  nuclear  test  site  at  Lop  Nor.  In 
this  project  we  have  been  able  to  construct  and  validate  attenuation  models  for  the  crust 
and  upper  mantle  in  western  China.  Previous  studies  have  demonstrated  the  great 
complexity  in  crustal  and  uppermost  mantle  seismic  velocity  structure  and  propagation 
characteristics  beneath  much  of  the  Qinghai-Tibet  plateau  and  the  surrounding  regions. 
Regional  seismic  phases  show  strong  lateral  variability  in  traveltime,  amplitude,  and 
frequency  content  due  to  this  very  complex  structure.  Due  to  a  lack  of  two-dimensional 
seismic  array  data,  the  three-dimensional  seismic  structure  and  propagation 
characteristics  of  the  crust,  lithosphere  and  upper  mantle  have  been  largely  unknown.  We 
have  helped  to  remedy  this  problem  using  data  from  relatively  new  broadband  stations  in 
the  region,  which  unlike  the  previous,  linear,  seismic  transects,  covered  most  of  the 
northeastern  portion  of  the  Qinghai-Tibet  plateau  in  a  semi-regularly  spaced  grid.  Using 
data  collected  by  Pi’s  on  the  NETS  (NorthEastem  Tibet  Seismic  experiment)  and 
ASCENT  (A  Seismic  Collaborative  Experiment  of  Northern  Tibet)  arrays  (5/07-6/09) 
which  consisted  of  more  than  110  temporary  broadband  stations  in  the  northeast  Tibetan 
Plateau,  Ordos  Plateau,  and  Qaidam  basin,  and  using  a  comprehensive  set  of  techniques, 
we  have  developed  new  seismic  velocity  and  attenuation  models.  Furthermore  we  have 
integrated  data  from  the  MIT  temporary  array  along  with  data  from  the  Namche  Barwa 
arrays.  These  new  data  sets  allow  the  application  of  new  techniques  designed  to  create 
robust  models  of  seismic  velocity  and  attenuation  in  which  we  reliably  estimate  the 
absolute  amplitude  of  the  velocity  and  Q  variations  across  western  China. 

We  have  been  able  to  leverage  existing  resources  to  significantly  extend  our  work 
beyond  what  we  proposed  to  deliver  originally.  We  have  extended  our  regional  phase 
attenuation  work  to  include  Northeastern  China.  Northeast  China  is  a  tectonically  active 
continental  craton  that  has  been  re-activated  by  Pacific  plate  subduction  and  India- 
Eurasia  collision.  We  have  extended  the  work  on  Sn  propagation  in  eastern  Tibet  to 
regional  phase  Sn  wave  propagation  across  northern  China  where  high  frequency  Sn  is 
more  prevalent  and  thus  allows  us  to  study  in  more  detail  high  frequency  Sn  propagation. 
The  temperature  and  the  corresponding  anelastic  variations  in  the  uppermost  mantle  are 
determined  by  calculating  Sn  Q.  We  have  collected  waveform  data  recorded  by  127 
stations  from  140  earthquakes  in  a  rectangular  region  from  30°N  to  60°N  in  latitude, 
100°E  to  145°E  in  longitude.  We  have  obtained  models  of  Sn  Q  for  northeast  China  using 
two  methods:  two  station  method  (TSM)  and  the  reverse  two-station  method  (RTM).  The 
former  method  measurements  are  contaminated  with  any  relative  changes  in  site  effect 
between  two  stations,  whereas  the  latter  are  not  and  are  likely  more  robust  in  an  absolute 
sense.  The  inversion  results  show  high  Q  values  in  the  Songliao  basin,  indicating 
relatively  low  temperature  and  low  viscosity  for  its  uppermost  mantle.  Checkerboard 
resolution  shows  good  resolution  at  a  scale  of  2°  x  2°  within  most  of  the  study  area. 
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2.  INTRODUCTION 


2.1  Tectonic  Setting 

Western  China  around  the  Lop  Nor  test  site  has  been  the  foeus  of  many  studies  related  to 
seismic  wave  propagation  and  velocity  structure  [e.g.  Li  et  ah,  2008;  Tilman  et  ah,  2003; 
Mitchell,  1997;  Rodgers  et  ah,  1997;  Rapine  et  al.  2003b;  McNamara  et  ah,  1997].  Here 
we  summarize  some  of  the  current  state-of-the-art  knowledge  of  the  seismic  structure  of 
the  Tibetan  plateau  and  its  relationship  to  tectonic  hypotheses. 


Southern  Tibet 
(Lhasa  Terrane) 


Northern  Tibet 
(Qiangtang  Terrane) 


Figure  1.  Rayleigh  wave  phase  dispersion  curves  from  northern  and  southern  Tibet. 
Diamonds  indicate  group  velocities  and  circles  indicate  phase  velocities.  Notice  the 
divergence  in  the  dispersion  curves  at  approximately  60  second  period.  This  is  most 
likely  a  result  of  anomalously  slow  velocities  in  the  uppermost  mantle  beneath 
northern  Tibet  (Rapine  et  al,  2003b). 


Body  wave  and  surface  wave  studies  have  indicated  a  seismically  slow  and  highly 
attenuating  uppermost  mantle  for  much  of  northern  Tibet  [Ni  and  Barazangi,  1983; 
McNamara  et  ah,  1997;  Rapine  et  al.,  2003a].  In  contrast,  in  southern  Tibet  and  the 
Qaidam  basin  (north  of  Kunlun  fault),  colder  mantle  has  been  implied  [e.g.  Ni  and 
Barazangi,  1983;  Tilmann  et  ah,  2003;  Galve  et  al.,  2002;  Li  et  al.,  2008]  (Figure  2). 
Tilmann  et  al.,  [2003]  successfully  imaged  the  downwelling  Indian  lithosphere  to  depths 
near  400  km  beneath  central  Tibet.  This  result  along  with  the  recent  P-wave  tomography 
images  of  the  Tibetan  plateau  indicate  that  the  Indian  continental  lithosphere  underlies 
southern  Tibet  and  an  estimated  minimum  of  1000  km  of  post-collisional  convergence 
between  India  and  Asia  has  occurred.  The  Indian  lower  crust  and  lithosphere  underthrust 
the  Himalayas  and  southern  Tibet  without  significant  internal  deformation  [e.g.  Kind  et 
al.,  2002;  Tilmann,  et  al.,  2003;  Kumar  et  al.,  2006].  In  southern  Tibet  a  variety  of 
seismic  observations  [Makovsky  and  Klemperer,  1999;  Yuan  et  ah,  1997;  Kind  et  al., 
2002;  Tilmann,  et  ah,  2003;  Huang  and  Zhao;  2006;  Li  et  ah,  2008]  and  electrical 
conductivity  [e.g.  Wei  et  al.,  2001]  all  point  to  a  weak  low-velocity  middle  crust 
underlain  by  a  strong  Indian  lithosphere  [Nelson  et  al.,  1996].  The  predominant  cause  of 
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this  low-viscosity  and  low-velocity  channel  in  the  mid-erust  is  wet  melting  of 
the  subducted  Indian  sediments  and  metasediments  [Nelson  et  al.,  1996;  Cotte  et  al., 
1999;  Rapine  et  al.,  2003b].  Therefore,  there  is  a  large  viscosity  contrast  vertically,  and 
deformation  is  confined  to  a  low-viscosity  layer  that  may  force  lower  erustal 
deformation  to  be  primarily  horizontal  [Bendick  and  Fleseh,  2007]. 

Another  effect  of  the  mid-erustal  low  veloeity  zone  is  the  high  attenuation  of  Lg  waves  in 
southern  Tibet  [McNamara  et  al.,  1996;  Reese  et  al.,  1999;  Xie,  2002a;  Phillips  et  al., 
2000,  2001,  2005).  In  eontrast,  the  erust  beneath  northern  Tibet  (Qiangtang  terrane)  is 
characterized  by  a  monotonically  increasing  seismie  velocity  with  about  8%  lower 
average  velocity  than  a  normal  continental  erust  [Rapine  et  al.,  2003a]  (Figure  2).  Lg 
attenuation  studies  of  northern  Tibet  indicate  low  Lg  attenuation  (Q  less  than  100  at  1 
Hz).  The  reason  for  such  low  Lg  attenuation  is  not  eompletely  understood,  but  it  seems 
related  to  high  temperature  and  a  partially  melted  erust  of  northern  Tibet  [Haeker  et  al., 
2000;  Fan  and  Lay,  2003]. 
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Figure  2.  A  simple  tectonie  map  for  the  Tibetan  Plateau.  The  major  features  shown  are 
ITS-  Indus  Tsangpo  Suture,  BNS-Bangong-Nujiang  Suture,  ATF-Altyn  Tagh  Fault,  KF- 
Kunlun  Fault. 

A  consequence  of  the  subduetion  of  Indian  continental  lithosphere  in  central  Tibet  is  that 
the  downwelling  of  lithospheric  material  inevitably  would  entrain  neighboring 
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A  maximum  seismogenie  thiekness  of  25  km  [Langin  et  al.,  2003]  is  observed  for 
northern  Tibet.  The  lithosphere  thiekness  is  about  140  km  according  to  S-P  converted 
phases  [Kumar,  et  al.,  2006];  however,  it  eould  be  thinner  in  the  interior  of  the  northern 
Plateau.  Rayleigh  wave  phase  veloeity  inversion  indicates  a  lithosphere  thickness  of  only 
about  120-140  km  thick. 
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asthenospheric  materials  along  which  would  be  induced  a  downward  flow.  This 
downwelling,  combined  with  tentative  evidence  for  southward-directed  subduction  of  the 
Qaidam  block  along  the  northern  plateau  margin  [Tapponnier  et  al.,  2001;  Kind  et  al., 
2002],  will  likewise  induce  a  downward  flow.  A  deficit  of  asthenosphere  must  be 
compensated  for  by  a  focused  upward  return-flow.  Such  an  upward  flow  would  provide 
an  explanation  for  a  low- velocity  body  imaged  by  Wittlinger  et  al.  [1996]  and  the  recent 
volcanism  just  south  of  the  Kunlun  Mountains,  and  a  mechanism  for  heating  the  crust  and 
gradual  erosion  of  the  remaining  Asian  lithosphere  beneath  northern  Tibet.  Tilmann  et 
al.’s  (2003)  P-wave  tomography  image  offers  only  a  glimpse  of  the  mantle  structure 
beneath  the  southern  Qiangtang  Terrane. 

The  Qaidam  Basin  (Figure  3)  is  a  low  relief  intermontane  basin  that  is  probably  underlain 
by  relatively  thick  lithospheric  mantle.  The  sediments  in  the  basin  probably  are  up  to  14 
km  thick  (Wang  and  Coward,  1990).  The  southern  margin  of  the  Qaidam  basin  is  marked 
by  the  North  Kunlun  fault  zone  which  is  the  major  thrust  front  along  the  southern  edge  of 
the  basin  where  there  appear  to  be  very  abrupt  changes  in  crustal  thickness  (~15  km)  (e.g. 
Zhu  and  Helmberger,  1998).  The  northeastern  edge  of  deformation  of  the  Tibet  plateau  is 
arguably  represented  by  the  fold-thrust  belts  of  the  Qilian  Shan  and  appears  to  mark  a 
transition  from  relatively  warm  and  seismically  slow  lithospheric  mantle  of  the 
Qiangtang  Terrane  to  colder  and  seismically  fast  lithospheric  mantle  beneath  the  Qaidam 
basin  and  possibly  the  Qilian  Shan  fold  and  thrust  belts. 

Dominating  the  northwestern  margin  of  both  the  Qaidam  Basin  and  the  Tibetan  Plateau 
proper  is  the  Altyn  Tagh  fault  (Figure  3).  In  scale  (1200  km  long)  and  geometry,  the 
Altyn  Tagh  suggests  a  truly  intracontinental  equivalent  to  the  great  transform  faults  like 
the  San  Andreas.  The  Altyn  Tagn  consists  of  many  strands,  not  all  of  which  are  currently 
active  (e.g.  Cowgill  et  al.,  2003).  Thrust  structures  of  the  Qaidam-Qilian  terranes  show 
little  evidence  of  being  distorted  by  the  Altyn  Tagh,  suggesting  that  the  former  are 
stronger  than  the  fault  zone  (Yin  et  al.,  2008).  Although  some  have  argued  that  the  Altyn 
Tagh  is  essentially  a  crust  fault,  detached  from  the  underlying  mantle  by  a  decollement 
(e.g.  Burchfiel  et  al,  1989),  recent  seismic  tomography  has  been  used  to  argue  that  it 
penetrates  to  about  140  km  (Wittlinger  et  al.,  1996).  Quaternary  basaltic  volcanism  near 
the  western  Altyn  Tagh  also  suggests  a  deeply  penetrating  structure  in  that  region,  while 
the  lack  of  such  volcanism  adjacent  to  the  Qaidam  terrane  could  be  construed  as  evidence 
for  the  lack  of  deep  penetration  along  its  northeastern  reaches  (e.g.  Yin  at  al.,  2008). 

All  of  these  terrains  sutures  are  near  the  Lop  Nor  test  site,  which  is  directly  to  the 
northwest  of  our  array.  This  results  in  an  incredibly  heterogeneous  lithospheric  seismic 
velocity  structure  that  represents  a  challenge  not  only  to  travel  time  calibration  and 
seismic  event  determination  but  also  to  amplitude  calibration  efforts.  In  this  project  we 
have  greatly  improved  the  ability  to  predict  both  travel  times  and  seismic  amplitudes  for 
the  purposes  of  seismic  event  discrimination  in  western  China. 
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Figure  3.  A  plot  of  seismic  stations  in  the  NETS(blue)  and  ASCENT  (green)  broadband 
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Figure  4.  A  plot  of  the  NECESS  broadband  array. 


Approved  for  public  release;  distribution  is  unlimited. 


5 


2.2  Attenuation  (0)  Tomography 

Attenuation  tomography  uses  estimates  of  path-averaged  attenuation  to  delineate  laterally 
varying  Q.  This  includes  most  spectral  methods  (Zor  et  ah,  2007).  An  example  of 
attenuation  tomography  for  Asia  using  coarsely  spaced  stations  is  shown  in  Figure  4  (Xie 
et  ah,  2006).  We  have  used  this  method  to  create  phase  blockage  maps  for  regions  of  the 
Middle  East  (Sandvol  et  ah,  2001). 
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Figure  5.  Lg  Q  estimates  using  spectral  methods  and  attenuation  (Q)  tomography  (Xie  et 
ah,  2006).  This  map  is  accompanied  by  a  map  of  frequency  dependence  (not  shown). 

Attenuation  tomography  relies  on  two-station  methods  for  measuring  inter-station  Q, 
which  cancel  the  source  effect  by  using  station  pairs  aligned  with  the  source.  These 
methods  are  relatively  straightforward  to  apply,  and  are  insensitive  to  source  radiation 
pattern  effects,  but  alignment  requirements  reduce  the  amount  of  data  that  can  be  used. 
We  measure  attenuation  with  the  two-station  (see  Xie  and  Mitchell,  1990;  Xie,  2002b; 
and  Zor  et  ah,  2007  for  descriptions  of  the  method,  including  rigorous  error  analysis)  and 
reversed  two  station  (Chun  et  ah,  1987)  methods.  The  reversed  method  eliminates  site  as 
well  as  source  effects.  Qo  and  r\  values  measured  using  these  methods  are  of  very  high 
quality  because  they  are  not  subject  to  the  trade-off  between  source  and  path  parameters, 
and  the  smaller  data  sets  can  be  easily  scrubbed  with  the  most  rigorous  quality  control. 
Measured  Qo  and  q  values  will  be  used  to  derive  a  long-wavelength  Q  map  for  the  study 
region.  This  model  can  be  subsequently  used  as  a  priori  knowledge  during  amplitude 
tomography  calculations.  Alternately,  inter-station  Q  measurements  can  be  included  as 
constraints  in  the  amplitude  tomography.  Figure  5  shows  a  preliminary  Lg  Qo  attenuation 
map  using  two  station  spectral  ratios.  Once  again,  ASCENT  is  ideally  designed  to 
tomographically  map  Q  in  this  region  with  a  high  degree  of  resolution,  and  will  increase 
resolution  in  surrounding  areas  when  included  in  continental  scale  studies. 
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Amplitude  tomography  uses  raw  amplitudes  rather  than  Q  estimates,  assumes  an  isotropic 
source  radiation  pattern,  and  solves  for  source  and  site  terms  in  addition  to  2-D 
attenuation.  An  example  of  amplitude  based  tomography  for  Asia  is  shown  in  Figure  6 
from  Phillips  et  al.  (2005).  The  amplitude  tomography  method  has  proven  successful  and 
performs  well  even  for  lower  quality  data  sets  such  as  the  surface-wave  amplitudes 
archived  by  the  International  Seismological  Centre  (Hearn  et  ah,  2008).  Amplitude 
tomography  does  not  require  any  special  station  geometry  and  can  utilize  many  more 
raypaths.  Constraints  from  high  quality  path  averaged  Q  such  as  those  described  above 
can  be  incorporated  into  the  inversion. 


Figure  6  Q  estimated  using  1-Hz  Lg  amplitudes  obtained  from  digital  waveform  data 
(Phillips  et  ah,  2005).  White  contours  represent  the  model  resolution  at  levels  of  0.1,  0.2, 
0.3,  and  0.4. 

Attenuation  causes  substantial  deviation  of  station  magnitude  from  the  true  event 
magnitude.  We  have  estimated  typical  ML  magnitude  corrections  using  our  past 
tomography  results  based  on  1  Hz  peak-energy  in  the  regional  waves.  These  maps  show 
that  regional  path  corrections  can  be  used  to  account  for  up  to  a  full  magnitude  unit. 
Station  gain  terms  can  cause  as  much  as  half  an  additional  magnitude  unit  correction. 
These  are  significant  corrections  from  a  monitoring  viewpoint.  Longer  period  waves, 
such  as  10-20  seconds  surface  waves,  produce  proportionately  smaller  corrections  due  to 
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the  influence  of  period  on  attenuation.  Thus,  the  10  second  surface  waves  used  to  image 
surface-wave  Q  across  China  by  Hearn  et  al.  (2008),  produce  corrections  of  only  a  tenth 
of  a  magnitude  unit  at  most,  although  the  station  corrections  can  again  cause  shifts  of  up 
to  half  a  magnitude  unit.  Long-period  amplitudes  are  less  sensitive  to  Q  variations  (Fisk, 
2006,  2007;  Fisk  et  al.,  2008)  than  short-period  phases.  Short-period  phases  however  are 
essential  to  regional  discrimination,  and  their  use  requires  accurate  attenuation 
corrections. 

3.  METHODOLOGY 

We  have  used  a  large  number  of  different  techniques  to  construct  attenuation  and  velocity 
images  of  eastern  Tibetan  and  northeastern  China;  therefore  it  is  necessary  to  discuss 
these  approaches  concisely. 

3.1  Two  Plane  Wave  Tomography  (TPWT) 

Seismic  surface  waves  travel  parallel  to  earth’s  surface,  sampling  different  depths  with  a 
range  of  frequencies.  Fundamental  mode  Rayleigh  and  Love  waves  are  most  sensitive  to 
earth  structure  at  depths  corresponding  to  ~  1/3  and  ~l/4  of  their  dominant  wavelengths, 
respectively  (Knopoff  et  al.,  1973).  Surface  waves  are  dispersive,  i.e.  relatively  longer 
period  surface  waves  have  greater  wavelengths,  allowing  them  to  travel  through  deeper 
parts  of  the  earth  at  greater  wave  speeds.  Using  a  range  of  frequencies,  we  can  take 
advantage  of  surface  wave  dispersion  to  calculate  phase  and  group  wave  speeds  as  a 
function  of  depth.  These  characteristics  of  surface  waves  make  them  valuable  tools 
for  studying  earth  structure  in  the  crust  and  uppermost  mantle. 

One  of  the  earliest  attempts  at  measuring  surface  wave  dispersion  was  the  single  station 
method.  This  method  uses  only  well-dispersed  records,  relying  heavily  on  the  initial 
phases  of  waves  generated  by  an  earthquake.  For  this  reason,  the  single  station  method 
requires  precise  knowledge  of  source  parameters.  Due  to  this  limitation,  it  is  often 
unreliable  for  measuring  surface  wave  dispersion.  A  more  traditional  and  widely  applied 
technique  is  the  two-station  method.  This  method  adapts  an  approach  similar  to  travel 
time  tomography,  employing  relative  travel  times  of  surface  waves  between  two  stations, 
lying  approximately  on  the  same  great  circle  path  (within  the  range  of  ±15°).  One 
advantage  of  this  approach  is  ease  of  implementation,  as  the  forward  problem  of  surface 
wave  propagation  is  reduced  to  seismic  ray  theory,  which  assumes  seismic  waves  consist 
of  infinite  frequencies.  However,  surface  waves  have  finite  frequencies.  Therefore, 
traditional  surface  wave  tomography  breaks  down  when  the  wavelengths  of 
heterogeneities  are  similar  to  wavelengths  of  surface  waves  (Li,  2011;  Zhou  and  Murphy, 
2005).  Additionally,  lateral  heterogeneities  between  events  and  seismic  arrays  may  cause 
scattering  and  multipathing.  Consequences  of  these  finite  frequency  effects  include 
distortion  of  ray  paths  and  variations  in  amplitudes  (Forsyth  and  Li,  2005).  For  these 
reasons,  amplitude  information  is  ignored  in  traditional  methods.  Moreover,  since  the 
traditional  method  requires  stations  to  be  on  the  same  great  circle  path  as  events,  not  all 
the  station-event  pairs  can  be  included  in  tomographic  inversions.  Therefore,  resolution 
and  robustness  of  phase  velocity  measurements  are  highly  affected  by  number  and  spatial 
distribution  of  ray  paths. 
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The  Two  Plane  Wave  Method  TPWT  (Forsyth  et  ah,  1998)  does  not  share  these 
limitations.  Rather,  it  uses  both  amplitude  and  phase  variations,  assuming  that  distortion 
of  a  wavefront  at  any  location  within  a  seismic  array  (i.e.  receiver  locations  and  inversion 
nodes)  can  be  expressed  as  the  sum  of  two  plane  waves  (Forsyth  et  ah,  1998;  Forsyth  and 
Li,  2005).  This  assumption  can  be  expressed  as  follows: 

)+  i^2  0xp(-i*=02),  (1) 

where  U  is  the  displacement  at  station  caused  by  event,  and  A  and  (p  are  amplitudes 
and  phases  of  two  plane  waves,  respectively.  Further,  the  phase  (0)  terms  in  Equation  (1) 
can  be  expressed  as: 

=  “^1,2  +  -  “t)  (2) 

where  represents  phases  of  first  and  second  plane  waves  at  a  reference  station,  and 

is  the  average  slowness,  ~  is  the  difference  between  travel  times  along  the  great 

circle  path,  measured  from  the  edge  of  study  area  to  the  A:*  and  reference  station. 
denotes  deviation  from  the  great  circle  path  for  both  plane  waves  (Forsyth  and  Li,  2005). 
Hence,  the  wavefield  at  each  station  for  a  single  frequency  can  be  described  by  three 
unknowns  for  each  plane  wave:  relative  amplitudes  and  phases  with  respect  to  a  reference 
station,  and  deviation  angles  from  the  great  circle  path. 


Figure  7.  A  schematic  of  local  Cartesian  coordinate  system.  The  reference  station  is  chosen 
using  the  minimum  RMS  misfit  in  amplitude  with  respect  to  mean  value,  x-axis  is  along  the  great 
circle  path  (gcp)  with  the  source,  andy-axis  is  90°  counterclockwise.  Triangles  show  stations. 

The  distance  between  the  source  and  the  center  of  the  seismic  array  is  greater  than  20  °. 
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To  invert  for  these  unknowns,  a  local  Cartesian  coordinate  system  is  set  up  for  each  event 
(Figure  7),  where  the  reference  station  is  the  origin,  the  x  direction  is  along  the  great 
circle  path  between  source  and  reference  station,  and  the  y  direction  is  90° 
counterclockwise  from  x.  The  reference  station  is  determined  using  the  lowest  root  mean 
square  error  of  amplitudes  with  respect  to  the  mean  value.  Alternate  methods  for  the 
reference  station  decision  include  using  the  shortest  event-receiver  distance  or  maximum 
amplitude. 
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Figure  8.  Map  views  of  amplitude  and  phase  sensitivity  kernels  for  20  mHz  Rayleigh 
waves,  obtained  using  the  Bom  approximation.  Panels  on  the  right  side  show  cross- 
sections  along  -1200  km  (white  solid  line  on  map  views).  The  hypothetical  station  is 
located  at  the  center,  and  source  is  to  the  west. 


Following  the  standard  procedure,  we  calculated  the  phase  wave  speeds  using  both 
Rayleigh  (vertical  component)  and/or  Love  (transverse  component)  waves  for  a  set  of 
frequencies  at  defined  points  within  a  study  area.  Then,  we  used  these  ID  dispersion 
results  for  obtaining  a  3D  shear  wave  stmcture.  More  specifically,  we  employed  the 
DISPER80  algorithm  (Saito,  1988)  for  each  individual  dispersion  curve  to  obtain  a  ID 
shear  wave  profile.  Next,  these  shear  wave  profiles  are  combined  to  constmct  a  3D 
model.  To  account  for  finite  frequency  effects,  we  incorporate  sensitivity  kernels  at  the 
phase  velocity  inversion  stage. 

In  order  to  account  for  the  finite  frequency  effects,  earlier  versions  of  TPWT 
implemented  Gaussian  kernels.  These  kernels  do  not  accurately  represent  finite  frequency 
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effects  as  they  assume  a  constant  Gaussian-shaped  sensitivity  function  perpendicular  to 
the  ray  path  (Forsyth  et  ah,  1998).  On  the  other  hand,  modified  versions  of  TPWT  use 
single  scattering  (Bom)  approximation  sensitivity  kernels  (Yang  and  Forsyth,  2006)  of 
Zhou  et  al.  (2004)  (Figure  8).  2D  sensitivity  kernels  for  fundamental  mode  surface  waves 
using  the  Bom  approximation  can  be  described  as  (Zhou  et  ah,  2004): 


K{T,oi)  = 


exp[— i[fe(A'  +  A"  —  A)  I2  —  n)  +  ^/4]} 


(3) 


where  K  is  the  sensitivity  kernel  as  a  function  of  propagation  direction  (r)  and  angular 
frequency  (co).  Single  and  double  primed  symbols  denote  waves  from  source  to  scatterer, 
and  from  scatterer  to  receiver,  respectively.  Symbols  with  no  primes  indicate  direct 
waves  from  source  to  station.  The  integer  n  is  the  polar  passage  index,  indicating  how 
many  times  a  wave  train  passes  thorough  the  source  or  its  antipode.  A  shows  the  distance 

between  source-scatterer  (^  ),  scatterer-receiver(^  ),  or  source-receiver(^)  .  Source  and 
receiver  polarization  terms  are  indicated  using  S  and  R,  respectively. 

The  real  and  imaginary  parts  of  Equation  (3)  give  amplitude  and  phase  kernels, 
respectively.  These  are: 

ai)  =  lm{K (r,  &>)}  ^4^ 

ft>)  =  -iie  {if  (r,  cj)}  (5^ 

In  the  above  equations,  scattered  and  direct  arriving  waves  are  due  to  the  same 

earthquake.  Therefore,  source  terms  (■^  ^  )  can  be  neglected.  In  our  case,  surface 

wave  tomography  is  a  minor  arc  problem.  Therefore,  a  wave  train  never  passes  through 
the  source  or  its  antipode  again.  Hence,  the  polar  passage  index  is  zero  (n  =  0),  and 
vanishes  from  Equation  (3).  Moreover,  in  the  case  of  Rayleigh  waves,  only  vertical 
components  of  seismograms  are  used.  For  this  reason,  receiver  polarization  vectors  {R) 
for  direct  and  scattered  waves  are  equal  (Yang  and  Forsyth,  2006).  However,  Love  wave 
energy  is  concentrated  on  the  transverse  plane;  and  polarization  terms  are  functions  of 
scattering  angle  (Snieder,  1986). 

In  reality,  sensitivity  kernels  are  not  the  same  for  the  whole  earth.  Rather,  they  are 
functions  of  earth  stmcture  and  may  change  along  the  propagation  path  in  their 
amplitudes  and  geometry.  However,  in  regional  tomography  studies,  research  areas  are 
relatively  small.  Most  probably,  variations  in  sensitivity  kernels  have  a  minor  role  in 
tomography  results,  and  therefore  1  make  a  final  assumption  and  neglect  these  variations. 
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3.2  Pn  and  Sn  Travel  Time  Tomography 

For  each  layer,  we  followed  the  Pn  tomography  technique  developed  by  Hearn  and  Ni 
(1994).  The  uppermost  layer  of  the  mantle  is  divided  into  a  set  of  two-dimensional  cells, 
in  which  the  slownesses  (inverse  velocities)  are  estimated.  The  Pn  travel  time  residuals 
are  described  approximately  by  the  perturbation  equation 

tij  —  Cli  “t  bj  “t  dijkSk 


where  ty  is  the  travel  time  residual  for  event  j  and  station  i,  a,  is  the  static  delay  for  station 
i,  bj  is  the  static  delay  for  event  j,  dyk  is  the  travel  distance  of  ray  ij  in  mantle  cell  k,  and  Sk 
is  the  slowness  perturbation  for  cell.  The  unknown  quantities  in  equation  (6)  are  the 
mantle  slowness  perturbation  Sk  and  the  station  and  event  delays  a,  and  bj,  respectively. 
The  station  and  event  delays  accommodate  variations  in  crustal  velocity  and  thickness. 

For  the  inversion,  we  used  a  cell  size  of  0.5°  x  0.5°.  The  slowness  values  in  each  cell 
were  resolved  using  the  LSQR  algorithm  (Paige  and  Saunders,  1982;  Yao  et  al.,  1999). 
While  we  attempted  to  resolve  azimuthal  anisotropy,  we  could  not  obtain  sufficient 
resolution  to  tightly  constrain  azimuthal  anisotropy  fast  directions  and  magnitude,  in  this 
project,  we  have  focused  on  the  Pn  and  Sn  velocity  variations  and  the  velocity  ratio, 
assuming  a  spatially  varying  isotropic  velocity  model.  It  is  important  to  note  that  one  of 
the  difficulties  in  measuring  Pn  and  Sn  azimuthal  anisotropy  is  the  difference  in  sampling 
depths  for  different  epicentral  distances.  Often  this  can  bias  these  types  of  measurements 
requiring  that  a  relatively  narrow  range  of  distances  be  used  to  solve  for  Pn  or  Sn 
anisotropy.  Typically,  including  the  anisotropic  terms  (cosine  and  sine  terms)  (Lii  et  al., 
2012),  has  a  marginal  impact  on  the  major  isotropic  velocity  variations.  Therefore,  our 
tomography  can  most  likely  provide  a  reliable  uppermost  mantle  velocity  structure  of  the 
study  area. 


3.3  Two  and  Reverse  Two  Station  0  Tomography 

The  amplitude  of  a  seismic  wave  may  be  described  by  an  exponential  attenuation 
equation  that  accounts  for  both  geometric  spreading  and  attenuation. 


4(/)  =  IXf)SsXf)Sj{f)G{^,j)txv 


4^ij 

v0(/) 


(V) 

where  A  is  the  observed  amplitude  between  source  (indices  j)  and  receiver  (indicesi)  for 
a  wave  of  frequency  f  recorded  at  distance  A.  Here,  f  is  the  instrument  response,  Sj  is 
the  source  amplitude,  Ssi  is  the  site  amplification  response,  and  v  is  the  wave  group 
speed.  Geometric  spreading  is  defined  as  the  following: 


G(A)  = 


1 


(8) 


For  cylindrical  spreading  m  is  0.5  and  for  spherical  spreading  m  is  1.0.  When  dispersion 
is  accounted  for,  these  can  increase  by  up  to  0.3  depending  on  signal  bandwidth.  We 
assume  this  function  is  frequency  independent;  however,  this  may  be  a  poor  assumption 
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in  many  cases  and  these  differences  will  be  mapped  into  our  effeetive  Q  or  more 
speeifically  the  seattering  component  of  our  effective  Q.  Spreading  relations  are 
important  but  poorly  documented  and  must  be  given  explieit  attention  in  attenuation 
problems.  The  attenuation  quality  factor,  Q,  can  be  assumed  to  be  frequeney  dependent, 
Q=Qo  f  with  Qo  being  the  attenuation  quality  faetor  at  1  Hz,  f  the  wave  frequency,  and  q 
the  frequency  dependence.  Values  for  Qo  and  q  depend  on  the  type  of  wave  used  but 
generally  q  lies  between  zero  and  one.  This  equation  does  not  aeeount  for  radiation 
pattern,  focusing,  or  anisotropic  effects.  However,  with  good  azimuthal  eoverage,  these 
effeets  should  average  out.  Both  Q  and  amplitude  tomography  methods  are  based  on  this 
equation,  with  appropriate  modifications. 

The  TSM  was  presented  by  Mitehell  (1995)  and  has  been  widely  used  (e.g.,  Nuttli,  1986). 
Xie  and  Mitehell  (1990)  use  the  assumption  of  spreading  sueh  as  that  shown  in  equation 
8  and  that  Q(f)=Qo^: 


(l  -  >7)111/- In  = 


V  A  A,(f)i,(f)d’: 

WFA)  U,{f)i.(f¥ 


/  (9) 

where  dj-di=A.  The  real  situation  for  applying  the  TSM  is  more  eomplieated  beeause  the 
perfeet  alignment  geometry  is  typically  not  obtainable,  especially  for  passive  seismic 
experiments.  In  practice,  event  to  station  paths  differ  by  a  small  angle  SO.  A  detailed 
analysis  of  this  angle  has  been  presented  by  Xie  et  al.  (2004).  Systematie  errors  could  be 
introduced  into  Qo  and  rj  values  determined  by  the  TSM  beeause  of  effects  of  attenuation 
out  of  the  path  and  anisotropic  source  radiation  patterns.  These  errors  ean  be  minimized  if 
a  threshold  value  dO^ax  is  used  to  limit  the  angle  50.  Xie  et  al.  (2004)  used  a  A^max  of 
±15°  in  their  study,  which  had  been  estimated  by  Der  et  al.  (1982).  Xie  et  al.  (2004)  also 
derived  an  equation  to  minimize  an  error  related  to  the  variation  of  inter-station  distanee 
at  1  Hz.  This  equation  (Xie  et  al.,  2004)  is  reorganized  here: 


Q  nf{d.-d,) 

where  dj-di  is  the  inter-station  distance,  SQo  is  the  error  in  the  measured  Qo  ,  and  5x  is  the 
total  unaeeounted  for  error  in  equation  (3),  which  assumes  an  ideal  geometry  and  ID 
structure.  In  this  study,  we  follow  Xie  et  al.  (2004)  by  setting  50.  to  15°  and  5x  smaller 
than  0.4.  Also,  the  equation  indicates  that  to  reduce  the  Q  error,  it  is  preferable  to 
maximize  the  inter-station  distance. 
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Figure  9  Schematic  drawing  of  the  geometry  of  RTS  and  RTE.  RTS-  (a)  the  ideal 
situation  where  the  source  a  (star  on  the  right),  the  station  i  and  j  (triangles),  and  source 
b  (star  on  the  left)  are  aligned  along  a  great  circle,  and  (b)  the  more  common  situation 
where  the  azimuth  difference  angles  from  sources  a  and  b  to  the  stations  are  denoted  by 
50a  and  56b  respectively.  RTE-  (c)  the  ideal  situation  where  the  station  i  (triangle  on  the 
right),  the  sources  a  and  b  (stars),  and  station  j  (triangle  on  the  left)  are  aligned  along  a 
great  circle,  and  (d)  the  more  common  situation  where  the  azimuth  difference  angles  from 
stations  i  and  j  to  the  events  are  denoted  by  50i  and  56,  respectively. 

The  main  disadvantage  of  the  TSM  is  caused  by  the  remaining  instrument  response  I  and 
neglected  site  response  Ss(0  The  instrument  responses  are  only  theoretically  identical  for 
stations  with  identical  seismometers  and  data  loggers.  In  practice,  instrument  calibration 
is  commonly  impossible  during  deployment,  which  also  does  not  preclude  the  influence 
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of  various  time-varying  uncalibrated  instabilities  on  the  recorded  amplitude.  It  is  also 
difficult  to  guarantee  the  operability  of  calibration  of  all  instruments.  The  site  effects 
have  been  suggested  to  have  strong  lateral  variations  and  are  likely  associated  with 
shallow  geological  structures,  especially  in  tectonically  active  zones  (Wald  and  Allen, 
2007;  Pasyanos  et  ah,  2009).  Therefore,  the  site  effects  independent  of  instruments  can 
contaminate  the  measurements  of  Q  if  we  use  the  TSM.  It  is  probably  an  over¬ 
simplification  to  assume  the  complicated  instrument  response  and  site  response  with  a 
merely  theoretically  consistent  model  of  instrument  response  and  a  neglected  site 
response. 


The  RTM  was  suggested  to  avoid  the  effects  of  neglected  Ss  terms  and  inaccurate  /  terms 
in  the  TSM  by  involving  one  more  event.  The  RTM  was  initially  suggested  by  Chun  et  al. 
(1987).  Figure  9  shows  the  geometry  of  the  RTM  including  its  two  cases:  the  Reverse 
Two  Station  (RTS)  paths  and  Reverse  Two  Event  (RTE)  paths.  The  ideal  case  for  RTS  is 
that  both  events  are  aligned  with  the  inter-station  path,  and  requires  the  four  epicentral 
distances  involved  to  be  within  a  regional  distance  range.  In  such  a  situation,  we  use  Aai, 
Aaj,  Abi,  and  Atj  to  denote  spectral  amplitudes  of  Lg  recorded  at  stations  i  and  j  for  events 
a  and  b,  and  dau  daj,  dbi,  and  dbj  the  corresponding  distances.  The  four  spectral  amplitudes 
can  be  expressed  as: 


AAf,dJ  =  SSf)RM,<P)Iiif)SsAf)GM^M 


4dai 

va_ 


A(f,dJ  =  S,{f)R,{f,(p)IXf)SsAf)GbAdbi)^^V 


^dgj 

¥dgi 


(11) 


Tzfdj^j 


If  Aai  is  divided  by  Aaj  and  Abi  is  divided  by  Abj,  we  get  (note  the  exponents  are  negative): 


A„..  R  7.  S,,  G„,. 


a  a  i 


Sg  Rg  Ij 

Sb  Rb  Ij  Ss,  G, 


exp 


^  4dgj  Tfd^ 


exp 


I  S,.  G 
7.  G„ 


-exp 


4dgj 

y^jQj 


4dgi 
^iQi  y 


dj  Ssj  G, 


exp 


Tfdbj  4db 


\ 

^  ■ (12) 


K^iQi 


v,a 


Like  the  TSM,  we  assume  that  the  velocity  structure  is  one-dimensional  and 
apparent  Q  values  are  identical  for  the  path  between  stations  i  and  j.  We  divide  the  two 
ratios  in  (1 1),  substitute,  and  obtain 


dgjdbj  ^ 

ydgjdi,.  J 


exp 


^{d. 

vQ^ 


'dbj+dbi) 


(13) 
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where  all  of  the  Q  parameters  above  are  a  function  of  frequency  and  m  is  the  exponent 
from  equation  8.  Therefore,  the  inter-station  (between  i  and  j)  apparent  HQ  value  can  be 
derived  as 


Q  4(daj-d,i-d,j+d,,) 


In 


AAj 


AjAi 


dqjdhj 
\dajdbi  j 


(14) 


whose  reciprocal  is  the  inter-station  apparent  Q.  Like  the  TSM,  (11)  to  (14)  are  all  in  the 
frequency  domain.  The  1/g  as  a  function  of  the  frequency /is 


1 


Qif) 


In 


MfK(f) 


d  du  ^ 

ai  bj 


(15) 


(15)  shows  that  the  RTS  does  not  require  any  assumptions  on  the  instrument  responses 
and  site  responses.  (15)  is  similar  to  (16),  but  it  requires  four  spectra  and  four  distances, 
not  the  two  spectra  and  two  distances  in  the  TSM.  A  pre-determination  of  m  is  also 
required,  as  it  is  with  the  TSM. 

Figure  9c  shows  that  the  geometry  of  ideal  RTF  closely  resembles  that  of  RTS  except 
that  the  positions  of  stations  and  events  are  switched.  This  case  was  not  presented  until 
the  study  of  Fan  and  Lay  (2003).  Similar  to  (1 1),  we  can  write 

A'dai 

4,  (/.  d., )  =  S,  (/K  (/,  p)/,  (/)S„  (/)G„  fc  y 

’fdgj 

A,  (/.  )  =  S,  (fK  (/.  v)l:  (/fc,  (/)G«  k,  > 

¥dbj 

A  (f.  d„  )  =  s,  (f)R,  (/,  <p)!,  (/)«,,  (f)0„  (d„ 

Similar  to  the  RTS,  we  assume  that  the  velocity  structure  is  one-dimensional  and 
apparent  Q  values  are  identical  at  events  a  and  b  in  the  RTE  case,  which  leads  from  (15) 
to  (17),  the  same  equation  as  in  the  RTS  case.  The  RTM  allows  the  event  locations  to 
deviate  from  the  inter-station  great  circle  by  two  small  angles  50a  and  5db  in  the  RTS,  and 
the  station  locations  to  deviate  from  the  inter-event  great  circle  by  two  small  angles 
50i  and  50j  in  the  RTE.  We  set  maximum  values  of  ±15°  (the  same  as  that  in  the  TSM) 
for  those  four  angles.  And  the  procedure  of  RTM  is  also  very  similar  to  that  of  TSM.  We 
follow  steps  as  used  in  the  TSM,  and  use  (15)  as  the  equation  to  calculate  Q  using  a  linear 
regression. 

The  main  disadvantage  of  RTS  is  the  typically  poor  path  coverage  because  of  its  rigorous 
geometrical  requirement.  Thus  its  feasibility  strongly  depends  on  the  seismic  network.  To 
satisfy  the  geometry  of  RTS,  continental  ray  paths  of  regional  distances,  inter-station  or 
inter-event  distance  not  less  than  150  km,  and  angle  difference  within  ±15°  are  required 
for  all  the  two  stations  and  two  events.  The  RTE  provides  a  dramatic  advantage  over  the 
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rigorous  requirements  of  the  RTS.  The  stations  are  not  required  to  be  deployed  in  the 
study  area,  and  can  be  deployed  around  a  seismic  zone  within  regional  distances,  which 
greatly  increases  the  practicability  and  accuracy  of  the  generalized  RTM.  Particularly,  the 
RTE  can  become  very  efficient  in  regions  with  relatively  active  seismicity  and  a 
relatively  difficult  natural  environment  for  dense  seismic  networks.  The  disadvantage  is 
that  its  resolution  is  reduced  by  probable  earthquake  mislocations  in  event  catalogs. 


The  second  application  of  the  RTM  is  that  we  can  use  it  to  solve  relative  site  responses. 
From  equation  10  we  obtain 


A  A 


I  s,. 

I  S I 


y^ajGbj  J 


exp 


Adai  ,  Adbj  nfdbi 


+  ■ 


V.g.  V.Q. 


(17) 


Using  equation  17  we  can  solve  for  the  logarithm  of  the  ratio  of  the  site  response  which 
gives  us  a  relationship  for  the  relative  site  response: 
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=  In  — + 


daA^bi-dai-dbj 
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And  can  be  further  transferred  to 


In  5., .  -  In  5., .  =  In  —  +  - 
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Tomographic  inversions  generally  use  least  squares  algorithms  such  as  LSQR  (Paige  and 
Saunders,  1982).  Constraints  on  the  problem  are  of  great  importance  and  can  be 
introduced  using  regularization  and  damping  techniques.  For  example,  spatial  variations 
in  attenuation  are  regularized  using  first  or  second  difference  smoothing  constraints. 
Source  models,  such  as  MDAC  (Walter  and  Taylor,  2002),  can  also  be  included  in 
tomographic  inversions  for  Q.  Moments  and  comer  frequencies  can  be  damped  to  known 
levels  for  special  events,  or  regularized  to  follow  a  best-fit  scaling  model.  Resolution  and 
covariance  are  quantified  using  impulse  responses,  checkerboard  tests,  and  matrix 
inversion  techniques.  Such  estimates  assume  that  the  model  equations  describe  the 
physics  perfectly,  which  we  know  is  not  completely  trae.  Dense  networks  such  as  the 
Iranian  combined  networks  have  allowed  us  to  investigate  what  are  very  often  unmodeled 
effects,  such  as  the  source  radiation,  focusing,  and  medium  anisotropy  effects  mentioned 
earlier. 
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Attenuation  tomography  uses  estimates  of  path-averaged  attenuation  to  delineate  laterally 
varying  Q.  This  includes  most  spectral  methods  (e.g.  Zor  et  ah,  2007).  An  example  of 
attenuation  tomography  for  Asia  using  coarsely  spaced  stations  can  be  found  in  Xie  et  ah, 
2006.  Our  attenuation  (Q)  tomography  relies  on  two-station  methods  for  measuring 
inter-station  Q,  which  cancel  the  source  effect  by  using  station  pairs  aligned  with  the 
source.  These  methods  are  relatively  straightforward  to  apply,  and  are  insensitive  to 
source  radiation  pattern  effects,  but  alignment  requirements  reduce  the  amount  of  data 
that  can  be  used.  We  have  measured  attenuation  with  the  two-station  (see  Xie  and 
Mitchell,  1990;  Xie,  2002b;  and  Zor  et  ah,  2007  for  descriptions  of  the  method,  including 
rigorous  error  analysis)  and  reversed  two  station  (Chun  et  ah,  1987)  methods. 

As  our  derivation  of  the  RTM  terms  has  shown,  the  reversed  method  eliminates  site  as 
well  as  source  effects.  Furthermore  we  have  been  able  to  use  the  reverse  two  station 
approach  to  isolate  the  site  and  calculate  the  relative  site  effect.  Thus,  after  correcting  for 
geometrical  spreading,  we  were  able  to  isolate  relative  site  terms  and  thereby  create  maps 
of  relative  site  terms  for  much  of  the  eastern  Tibetan  plateau.  We  can  use  these  relations 
to  determine  both  path  based  Q  values  that  are  independent  of  the  site  effect  as  well  as  a 
relative  site  effect  term.  We  have  found  this  to  improve  attenuation  models  in  particular 
for  Pg  attenuation  in  eastern  Tibet.  This  has  been  previously  difficult  in  the  Tibetan 
plateau  because  of  a  lack  of  contemporaneously  recording  stations. 

Qo  and  r\  values  measured  using  these  methods  are  of  very  high  quality  because  they  are 
not  subject  to  the  trade-off  between  source  and  path  parameters,  and  the  smaller  data  sets 
can  be  easily  scrubbed  with  the  most  rigorous  quality  control.  Measured  Qo  and  r|  values 
are  used  to  derive  a  long-wavelength  Q  map  for  the  eastern  Tibetan  Plateau.  This  model 
can  be  subsequently  used  as  a  priori  knowledge  during  amplitude  tomography 
calculations.  Alternately,  inter-station  Q  measurements  can  be  included  as  constraints  in 
the  amplitude  tomography.  Once  again,  the  combined  ASCENT,  NETS,  and  Namche 
Barwe  arrays  are  ideally  designed  to  tomographically  map  Q  in  this  region  with  a  high 
degree  of  resolution,  and  have  increased  resolution  in  surrounding  areas  when  included  in 
continental  scale  studies. 

3.4  0  Anisotropy  Method 

Given  the  reliability  of  our  reverse  two  station  Q  estimates,  we  determined  that  we  should 
be  able,  given  sufficient  ray  path  coverage,  to  resolve  azimuthal  variations  in  Q.  We 
suggest  that  the  azimuthal  anisotropy  of  l/Qig  can  be  described  as 

l/Qig  =  A  +  Bs\n^{(p-0)+Cs\n^{(p-0),  (20) 

where  A  is  the  isotropic  coefficient  and  B  and  C  are  the  2nd  order  and  4th  order 
anisotropic  coefficients  respectively.  The  angles  (p  and  6  denote  the  path  azimuth  and 
high-g  direction,  respectively.  Furthermore,  (20)  can  be  transformed  to 
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(21) 
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We  notice  that  (2)  is  similar  to  the  expression  of  azimuthal  anisotropy  for  surface  waves 
(Smith  and  Dahlen,  1973).  Typically  the  4(p  terms  are  considered  negligible  (Montagner 
and  Nataf,  1986)  in  many  studies  (e.g.,  Hearn,  1996;  Huang  et  al.,  2003;  Yao  et  al., 
2010).  Moreover,  a  comparison  between  velocity  and  attenuation  azimuthal  anisotropy 
(Chapman,  2009)  and  our  observations  suggest  that  the  effect  from  4^  terms  can  be 
neglected  for  Qig  measured  from  the  vertical  component  signal.  Thus  (21)  will  be 
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where  A  is  the  isotropic  coefficient  (or  isotropic  HQig)  and  B  the  anisotropic  coefficient 
(or  anisotropic  magnitude). 


We  can  assume  that  A,  B,  and  6  are  all  frequency  dependent.  The  lateral  variations  in  A 
and  B  can  be  solved  from  4-spectral-amplitude  groups  at  different  frequencies  by 
dividing  the  study  area  into  a  two  dimensional  set  of  square  cells.  All  Lg  windows  are 
manually  picked  and  show  a  nearly  constant  vig  of  3.5  km/s,  and  the  study  of  Chapman 
(2009)  suggests  that  the  lateral  variation  in  velocity  is  probably  negligible  using  this 
method.  The  geometrical  spreading  factor  m  of  Lg  is  typically  set  at  0.5  (Hasegawa, 
1985).  We  assume  that  A,  B,  and  6  are  constant  inside  each  cell  at  each  frequency. 
Referring  to  the  method  of  Q  tomography  of  Xie  and  Mitchell  (1990),  for  each  4- 
spectral-amplitude  group  and  each  frequency/,,  (22)  can  be  finally  derived  as 
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^kn  =  An  ^\^kn  ^  4,  =  -^cos6>,„ ,  =  -^sin6>,„, 

where  dk  and  (p^  are  the  distance  and  azimuth  travelled  by  path  ij  in  cell  k,  respectively, 
and  Akn,  Bk„,  and  6^y  the  isotropic  coefficient,  anisotropic  coefficient,  and  high-g 
direction  of  cell  k  and  frequency/,,  respectively.  Akn,  Bkn,  and  6kv  are  unknown  variables 
and  are  assumed  to  be  frequency  dependent. 


The  linear  system  of  (23)  can  be  written  in  matrix  form  D=Gm.  We  denote  the  number  of 
4-spectra-amplitude  groups  by  Na  and  the  number  of  discrete  frequencies  by  A/.  The 
vector  D  has  its  elements  calculated  from  the  terms  on  the  right  side  of  the  equal  sign  in 
(4),  and  the  length  of  it  is  The  vector  m  consists  of  all  unknown  parameters 

including  A/,  Ykn,  and  Zkn  in  (23)  of  all  cells  and  frequencies,  and  the  length  of  it  is 
3MNNf  if  the  study  area  is  divided  into  M^N  cells.  The  matrix  G  is  a  large  and  sparse 
matrix  with  the  size  of  NaNf  ^  SMNNp-  We  use  an  LSQR  algorithm  (Paige  and  Saunders, 


Approved  for  public  release;  distribution  is  unlimited. 


19 


1982)  to  solve  this  inverse  problem.  Like  the  Pn  velocity  anisotropy  (Hearn,  1996)  and 
surface  wave  azimuthal  anisotropy  (e.g.,  Yao  et  al,  2010),  the  corresponding  elements  in 
the  solution  vector  m  can  be  used  to  determine  Ak„  and  Bk„,  and  the  high-g  direction  0,^. 
by: 


=  z 

B,.  =  2^r,l+zl 


2 

kn 


^kn~ 


\ 


1  c 

9  (3+-arcta4'" 

2  Wj 


1  8‘0+-arctb^'" 
2 


1 


—a  r  c  ta-« 


Z  ^ 

^kn 


\^kn  J 


if  Ykn^^ 
if  };„<0,  Z,„>0 


(24) 


\^kn  J 


if  7,„<0,  Z,„<0 


The  result  is  regularized  by  using  a  nine-point  spatial  smoothing,  which  was  used  in  Qig 
tomography  by  Suetsugu  and  Nakanishi  (1985)  and  Xie  and  Mitchell  (1990).  If  k(p,  q) 
denotes  the  ^h  cell,  where  p  and  q  are  cell  numbers  in  the  latitude  and  longitude 
directions,  respectively,  a  set  of  regularization  equations  are  organized  by  counting  all 
neighboring  cells  of  the  ^h  cell  to  build  a  matrix  L: 
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Two  damping  factors  A,\  and  Z2  are  applied  to  balance  the  trade-off  between  isotropic  and 
anisotropic  terms,  while  X\  is  for  Xkn  and  X2  is  for  and  Z^,.  Two  damping  factors  X\ 
and  X2  are  applied  to  balance  the  trade-off  between  isotropic  and  anisotropic  terms,  while 
X\  is  for  Xkn  and  X2  is  for  Ykn  and  Z^. 


We  have  also  continued  work  on  creating  a  reliable  tomographic  model  for  the  Pn 
velocity  structure  in  western  China  including  station  corrections  or  crustal  delay  times. 
Data  from  recent  INDEPTH  IV  and  other  broadband  deployments  in  Tibet  are  used  to 
invert  for  crustal  delays  and  Pn  velocities  beneath  the  eastern  Tibetan  Plateau  and 
surrounding  regions.  The  average  Pn  velocity  for  the  region  is  8.1km/s  but  varies  from 
7.8  to  over  8.3km/s.  Generally  low  Pn  velocities  are  found  in  northeastern  Tibet  in  the 
Qiangtang  and  Songpan-Ganzi  terrains.  This  includes  a  zone  of  very  low  Pn  velocity 
northwest  of  the  Longmen  Shan  thrust,  along  the  eastern  Kunlun  fault,  and  beneath  the 
eastern  Qilian  Shan.  A  region  of  high  Pn  velocity  underlies  the  eastern  end  of  the 
Bangong-Nujiang  Suture  (bounded  by  90“-98“E).  It  could  represent  part  of  the 

underthrusted  Indian  shield.  The  region  between  the  Qilian  Shan  and  Kunlun  Shan  is  also 
characterized  by  high  Pn  velocity  with  several  zones  of  extremely  high  velocity.  This 
includes  two  high  velocity  features  beneath  Qaidam  basin  and  Gonghe  basin.  These 
features  may  correspond  to  cratonic  fragments  that  accreted  during  the  closure  of  the 
Tethys  Ocean  and  have  impeded,  but  not  stopped,  the  northward  growth  of  the  plateau. 
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The  station  delays  show  the  thickest  crust  is  beneath  Tanggula  Shan  in  central  Tibet  and 
should  be  the  result  of  the  internal  deformation  within  Qiangtang  terrain  due  to  the 
collision  of  Indian  and  Eurasia  plates.  There  is  a  significant  decrease  of  crustal  thickness 
in  the  northeast  plateau  and  lateral  variation  also  exists  within  the  region. 

To  investigate  the  lateral  resolution  of  the  inversion  parameters,  checkerboard  tests  are 
used  to  determine  the  most  appropriate  damping  factors  and  the  approximate  size  of 
anomalies  which  can  be  retrieved.  Two  reference  values  for  ^  and  B  are  both  set  at  1/300, 
respectively,  and  the  perturbation  of  A  is  50%,  which  is  selected  to  approximately 
simulate  the  observation  of  anisotropic  attenuation  from  our  data  (Figure  10).  Two 
sinusoidal  functions  designed  for  north-south  (N-S)  and  west-east  (E-W)  synthetic  high  Q 
directions,  respectively,  are  used  to  calculate  the  variations  of  A  and  B.  Thus  each 
checkerboard  block  has  the  maximum  (1/200)  or  minimum  (1/600)  of  ^  and  maximum  of 
B  (1/300)  at  its  center,  and  the  average  of  A  (1/300)  and  minimum  of  B  (0)  along  all  the 
edges  of  the  synthetic  Q  anomalies.  Staggered  0  values  of  0°  and  90°  are  used  in 
neighboring  blocks.  The  checkerboard  tests  of  1°  x  1°  anomalies  are  likely  irretrievable, 
especially  for  the  anisotropic  terms,  which  is  probably  because  of  insufficient  path 
density.  We  further  test  four  checkerboard  patterns  where  all  designed  anomalies  have  the 
same  size  of  2°  x  2°  (Fig.  3).  Among  the  four  patterns,  Pattern  1  and  Pattern  2  have  the 
same  variations  of  A  and  B  but  opposite  6  distribution;  and  Pattern  3  and  Pattern  4  have 
the  opposite  variations  of  A  and  B  to  Pattern  1  and  Pattern  2,  and  the  same  0  distribution 
as  Pattern  1  and  Pattern  2,  respectively.  A  normally  distributed  random  noise  of  0%,  15%, 
and  25%,  is  added  into  all  the  tests.  We  simultaneously  solve  A,  B,  and  ^by  dividing  the 
region  into  0.5°  x  0.5°  square  cells  in  tomography.  The  ideal  damping  factors  Ti  and  A2 
should  be  selected  to  minimize  the  trade-off  between  lateral  variations  of  A  and  B,  which 
occurs  when  a  variable  Rja,  defined  from  the  residuals  of  A  and  B  as 


R. 


(26) 


is  equal  to  1,  where  Aq  and  Bq  denote  parameters  in  the  given  model  and  A  and  B  denote 
parameters  in  the  solution.  The  Rja  in  (26)  only  includes  residuals  in  the  area  defined  by 
limiting  the  hit-count  to  more  than  5  within  each  cell  of  the  model. 

An  example  of  how  the  best  Ti  and  A2  are  determined  for  Pattern- 1  in  Figure  10  is  shown 
in  Figure  11,  where  a  grid  search  is  used  for  different  A\  and  A2  pairs  varied  from  0.02  to 
0.3,  and  contours  of  Ria  value  are  illustrated  and  labeled.  An  important  observation  is  that 
an  increase  in  leads  to  an  increase  in  the  residual  of  A  and  a  corresponding  decrease  in 
the  residual  of  B,  and  an  increase  /I2  causes  an  increase  in  the  residual  of  B  and  a  decrease 
in  the  residual  of  A  low  residual  of  A  probably  deserves  priority  in  our  problem 
because  A  will  be  very  important  in  a  further  discussion  on  its  frequency  dependency. 
Thus  we  determine  the  best  Ai  and  A2  pairs  by  satisfying  a  sufficiently  low  residual  of  A 
and  a  point  within  a  15%  perturbation  of  Rja  of  1  (between  0.85  and  1.15)  among  the 
contours  in  Fig.  4.  Moreover,  we  observe  that  the  level  of  random  noise  may  slightly 
affect  Ria,  where  the  contours  change  substantially  by  increasing  the  noise  from  0%  to 
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25%,  but  this  is  not  as  significant  as  the  influence  from  the  ehanges  in  damping 
parameters.  We  also  observe  that  we  are  able  to  maintain  a  stable  and  relatively  high 
resolution  with  a  level  of  random  noise  of  up  to  15%,  but  a  higher  level  of  random  noise 
rapidly  deereases  the  resolvability  (defined  by  Zelt  [1998])  and  makes  the  inversion 
unstable.  We  finally  determine  that  X\  and  X2  can  be  typieally  seleeted  at  0.15  and  0.25, 
respectively,  based  on  a  number  of  sueh  tests  using  a  grid  search.  This  actually  means  that 
the  trade-off  between  A  and  B  is  eonditionally  minimized  to  less  than  15%.  The  results  of 
the  four  tests  indieate  good  resolution  for  anomalies  at  2°  x  2°  scale  with  15%  random 
noise  when  we  use  a  cell  size  of  0.5°  x  0.5°  (Fig.  3). 


Pattern- 1 


Pattern- 2 


Pattern-3 


Pattern-4 


92  96*  100’  92  96*  100’  92  96*  100’  92  96*  100’ 

0.002  0.003  0.004  0.005  “ 

Figure  10.  A  resolution  test  to  determine  whether  we  are  able  to  resolve  azimuthal 
anisotropy  given  our  reverse  two  station  path  eoverage. 
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Figure  11.  An  example  of  the  grid  search  result  on  the  co-influence  of  damping  factors  A.i 
and  A2  on  the  residual  ratios  of  isotropic  and  anisotropic  coefficients  A  and  B  when  (a)  no 
random  noise,  (b)  15%  random  noise,  and  (c)  25%  random  noise  are  added,  respectively. 
The  checkerboard  pattern  of  this  example  is  shown  in  Pattern  1  of  Figure  10.  The  gray 
color  backgrounds  represent  scales  of  the  residuals  of  A.  The  labeled  contours  represent 
the  distribution  of  Ria  defined  in  (26).  The  effect  from  added  random  noise  level  is 
dominant  over  that  from  damping  factors.  The  best  damping  factor  pair  is  selected  when 
the  residual  of  A  is  at  a  minimum  value  and  the  perturbation  of  Ria  is  between  ±15 
percent  (0.85  to  1.15).  This  shows  that  the  trade-off  between  isotropy  and  anisotropy  is 
conditionally  minimized  at  a  secondary  bias  for  the  approach  of  a  realistic  isotropic 
coefficient.  Typically,  Aj  and  A2  can  be  selected  at  0.15  and  0.25,  respectively,  based  on 
many  such  grid  search  results. 


3.5  Catalog  Amplitude  Tomography  (ML  tomography) 

In  general  we  have  taken  the  approach  of  using  multiple  independent  approaches  to 
measure  seismic  velocities  and  attenuation.  Thus  in  an  effort  to  test  or  validate  our 
crustal  attenuation  from  temporary  broadband  stations  in  the  region  we  have  begun  to 
mine  amplitude  data  from  the  official  Chinese  Earthquake  Bulletin.  The  data  from  the 
Chinese  Bulletin  was  also  used  to  estimate  crustal  attenuation.  This  was  done  in  a 
manner  similar  to  Hearn  et  ah,  2004.  That  study  used  data  from  the  Annual  China 
Bulletin  of  Earthquakes  but  was  limited  to  only  about  100  stations.  In  2008  China  began 
incorporating  their  sub-network  data  into  their  bulletin  and  now  has  over  1000  stations 
distributed  across  China.  Bulletin  data  lists  event  locations,  arrival  times  for  both 
regional  and  teleseismic  phases,  and  amplitude  and  period  data  for  the  purpose  of 
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estimating  ML,  MS,  MB,  and  Mb  magnitudes.  In  this  study  we  have  used  the  ML 
amplitudes  and  periods. 

The  China  Bulletin  is  an  amalgamation  of  national  network  data  and  28  sub-networks. 
Data  was  downloaded  from  the  online  website  of  China  Earthquake  Data  Center 
(http://data.earthquake.cn/data)  and  contains  data  for  events  as  small  as  magnitude  1.0. 
The  online  information  does  not  list  station  locations.  We  obtained  these  from  other 
sources  both  official  and  non-official.  We  also  mass  relocated  stations  using  the  source- 
receiver  distances  and  azimuths  in  the  bulletin.  Best  locations  were  chosen  on  the  basis 
of  the  residuals  and  verified,  when  possible,  by  the  use  of  Google  Earth.  Elevations  were 
found  using  an  API  (Application  Programming  Interface)  for  Google  Earth.  Because  the 
bulletin  is  a  combined  bulletin,  there  were  many  inconsistencies  and  blunders  in  the  data 
set  that  had  to  be  rectified. 

We  found  most  of  the  raw  data  to  be  useful  but  listed  source-receiver  distances  and 
azimuths  often  were  off.  This  occurs  because  the  sub-network  data  was  referenced  to  the 
sub-net  location  rather  than  the  national  location  of  the  bulletin.  In  one  case,  for  the 
Tianshan  sub-network,  the  azimuth  data  made  no  sense  at  all.  In  addition,  since  sub¬ 
networks  share  data,  there  are  sometimes  multiple  readings  for  each  station  and  it  is  not 
clear  which  was  used  for  the  final  location.  For  amplitude  phases  both  the  peak-to-peak 
amplitude  and  period  were  used  (the  reported  period  is  the  pulse-width  at  maximum 
amplitude).  Amplitudes  were  measured  using  both  velocity  and  displacement 
seismographs  and  designated  as  such.  However,  amplitudes  were  often  reported  using 
different  units  of  nanometers,  micrometers,  and  deca-nanometers  (10’^  meters)  without 
specifying  which  unit  it  was.  Sub-network  seismogram  data  was  reported  using  deca- 
nanometers  for  the  velocity  seismograph  readings  and  nanometers  for  the  displacement 
seismograph  readings.  National  network  data  was  reported  using  nanometers  through 
2010  and  micrometers  from  July  21,  2012  to  present.  Between  these  dates  the  units 
alternated  and  these  data  were  excluded.  We  adjusted  all  data  to  nanometers  before  data 
analysis;  however,  since  the  inversion  contains  station  gains  any  mistake  about  units  will 
be  absorbed  into  those  terms. 

3.6  Joint  Inversion  of  Receiver  Functions  and  Surface  Waves 

We  used  receiver  functions  from  seventy-eight  stations  from  the  ASCENT  and  NETS 
arrays  in  joint  inversions  for  the  velocity  models.  We  used  the  method  described  by  Julia 
et  ah,  2000  with  some  minor  modifications.  Both  phase  and  group  velocities  across  the 
eastern  plateau  were  used  to  interpolate  velocities  at  the  station’s  coordinates  using  a 
least  mean  squares  approach.  The  phase  velocity  model  was  taken  from  Ceylan  et  ah, 
2012  and  the  group  velocity  measurements  were  taken  from  Pasyanos,  2010. 

Receiver  functions  for  each  station  were  cut  from  5  seconds  before  the  P-wave  arrival  to 
65  seconds  after  the  first  arrival.  We  used  the  method  of  Liggoria  and  Ammon,  1999  to 
compute  the  time  domain  receiver  functions.  The  receiver  files  were  then  grouped  and 
stacked  according  to  their  ten  degree  distance  bins.  This  allowed  us  to  model  the  normal 
move-out  of  the  PS  Moho  phases  to  help  constrain  the  velocity  inversion.  The  events 
from  65°  to  85°  were  grouped  as  one  due  to  a  lower  number  of  files  at  these  distances. 
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The  events  from  35°  to  55°  were  put  into  5°  groupings  as  these  files  were  much  more 
numerous  at  these  distances  and  to  allow  for  easier  processing.  We  used  a  Gaussian  filter 
alpha  of  2.5  for  all  receiver  functions. 

We  chose  to  allow  for  more  receiver  file  influence  or  weighting  than  dispersion  (receiver 
function  =  0.6;  dispersion  data  =  0.4).  The  optimal  smoothing  parameter  was  found  to  be 
0.7,  which  gave  the  best  balance  between  producing  a  smoother  velocity  model  and 
allowing  for  a  reasonable  fit  to  the  receiver  functions.  We  compiled  all  receiver 
functions  and  interpolated  between  all  70  stations  in  order  to  produce  a  three  dimensional 
velocity  model.  The  velocity  models  were  then  used  to  construct  a  velocity  surface  at  50 
km,  80  km  and  100  km  to  spatially  evaluate  the  models  in  accordance  with  accepted 
theory. 

4.  RESULTS  AND  DISCUSSION 

4.1  Two  Plane  Wave  Rayleigh  Wave  Tomography 

We  used  fundamental  mode  Rayleigh  waves  and  two-plane  wave  tomography  (TPWT) 
for  phase  velocity  inversions  (Forsyth  and  Li,  2005).  We  chose  earthquakes  with 
epicentral  distances  between  20°  and  120°,  depth  <  100  km,  and  MS  >  5.7.  We  only  used 
the  vertical  component  in  order  to  avoid  shear  and  Love  wave  interference,  and  long 
period  noise  that  may  exist  on  horizontal  components.  The  events  with  low  signal  to 
noise  (S/N)  ratios  were  eliminated  visually.  Our  data  was  recorded  at  the  INDEPTH-IV 
array  (74  broadband  stations),  which  was  deployed  across  Northern  Tibet.  We  also 
included  data  from  the  Namche-Barwa  seismic  array  (Sol  et  ah,  2007)  in  southeastern 
Tibet.  Prior  to  data  processing,  all  instrument  responses  were  converted  to  the  first 
generation  STS2.  We  measured  phase  velocities  for  13  frequency  bands  ranging  between 
20  and  143  s,  and  utilized  10-mHz-wide  Butterworth  filters  centered  on  each  frequency. 
In  order  to  avoid  dependence  on  earthquake  magnitude,  we  normalized  the  seismic  wave 
amplitude  for  each  event,  then  measured  phases  and  amplitudes  of  visually  windowed 
data  using  Fourier  analysis.  We  used  a  total  of  174  earthquakes.  In  order  to  increase  the 
reliability  of  our  inversions,  we  only  included  events  that  were  recorded  by  at  least  10 
stations. 
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Figure  12.  The  difference  between  isotropic  phase  velocities  obtained  using  ambient 
noise  tomography  (ANT)  [Yang  et  ah,  2010]  and  the  two-plane  wave  approach  (TWPT) 
for  periods  of  a)  25,  b)  30,  and  c)  40  s.  The  third  column  shows  the  difference  between  the 
two  methods  in  km/s.  These  results  show  that  within  the  seismic  array  (with  the  exception 
of  Qaidam  Basin),  phase  velocities  obtained  by  two  methods  are  similar  (~±0.06  km/s). 
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Figure  13.  Ray  path  coverage  for  each  period  (T)  with  the  number  of  events  (N). 
Isotropic  phase  velocities  are  stable  for  T>=25  s  (the  periods  with  N>=60  earthquakes). 
White  triangles  show  the  seismic  station  locations.  For  periods  less  than  40  seconds, 
there  is  not  enough  coverage  on  northwestern  edge  of  our  array. 

The  first  step  in  our  inversions  is  calculating  isotropic  phase  velocities  for  each 
frequency.  Our  experiments  with  TPWT  show  that  ~60  events  are  enough  to  obtain  stable 
results.  In  order  to  better  constrain  our  inversion  parameters  (i.e.  smoothing  length  and 
damping),  we  compared  our  results  with  those  of  Ambient  Noise  Tomography  (ANT) 
(Yang  et  al.,  2010).  Although  both  methods  measure  surface  wave  dispersion,  ANT  and 
TPWT  are  fundamentally  different  and  use  independent  data  sets:  TPWT  implements 
surface  waves  created  by  earthquakes  while  ANT  excludes  earthquakes  and  uses  ambient 
noise  of  the  earth.  Additionally,  while  ANT  uses  ray  theory  to  obtain  phase  maps,  TPWT 
assumes  that  distortion  of  wavefronts  at  each  station  can  be  expressed  as  the  sum  of  two 
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plane  waves  (Forsyth  and  Li,  2005).  Due  to  the  short  wavelengths  of  ambient  surface 
waves,  ANT  can  typically  obtain  more  reliable  results  at  crustal  depths  than  TPWT. 
There  is  an  overlapping  frequency  range  between  these  two  methods  (in  our  case  20-50 
s). 

The  comparison  of  isotropic  phase  velocities  indicates  that  within  the  seismic  array,  these 
two  methods  reveal  very  similar  results  (differences  less  than  ±  0.06  km/s  with  the 
exception  of  Qaidam  Basin)  for  periods  between  25  and  40  s  (Figure  12).  With  the 
exception  of  125  s  (~170  km),  the  checkerboard  tests  (Figure  14)  imply  that  resolution 
within  the  seismic  array  is  ~  150  km  or  better  at  all  periods,  sufficient  to  image  all  of  the 
anomalies  discussed  in  this  study.  Due  to  insufficient  ray  paths  from  the  W-NW,  the 
checkerboard  tests  exhibit  evidence  for  some  lateral  smearing. 
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Figure  14.  The  resolution  test  results  for  different  periods  using  checkerboard  synthetic 
model.  The  anomalies  in  the  input  model  are  1.5  °x  1.5°.  The  periods  are  indicated  in  the 
upper-right  corner  of  each  panel. 
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Figure  15.  Anisotropic  phase  velocity  maps  for  12  periods  (22-143  s).  The  solid  black 
lines  indicate  azimuthal  fast  direction  and  amplitude  of  anisotropy.  Note  that  color  scales 
are  different  for  each  map.  The  regions  with  errors  of  2  cr  <  0.06  km/s  (where  a  is  the 
standard  deviation)  are  shaded.  HVBl  and  HVB2  show  the  high  velocity  bodies 
addressed  in  the  text. 
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Figure  16.  Shear  wave  velocities  (km/s)  for  depths  between  30  and  210  km.  Depths  are 
indicated  on  each  panel.  Note  that  each  panel  has  a  different  color  scale.  The  dark  gray 
lines  show  major  tectonic  features.  HVBl  and  HVB2  show  the  high  velocity  bodies 
addressed  in  the  text.  The  regions  with  errors  of  2(7  < 0.06  km/s  (where  a  is  the  standard 
deviation)  are  shaded,  based  on  peak  depth  sensitivity  of  corresponding  frequencies. 


The  second  step  of  phase  velocity  inversion  includes  determining  anisotropic  phase 
velocities  and  seismic  fast  directions  (Figure  15).  Unlike  traditional  ray  theory,  TPWT 
accounts  for  phase  and  amplitude  changes.  A  disadvantage  of  this  method  is  its 
sensitivity  to  finite  frequency  effects  such  as  scattering  and  wave  front  healing  (Hung  et 
al.,  2001;  Nolet  and  Dahlen,  2000).  To  overcome  these  limitations,  we  employed  2D 
sensitivity  kernels  [Yang  and  Forsyth,  2006]  for  both  amplitude  and  phase  perturbation 
for  each  period.  Azimuthal  anisotropy  was  simultaneously  solved  with  Rayleigh  wave 
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phase  velocities.  Because  the  shortest  wavelength  for  our  data  is  ~65  km  (for  20  s  period 
and  taking  the  average  phase  velocity  as  3.18  km/s),  we  used  a  constant  smoothing  length 
of  80  km  for  all  periods. 

In  order  to  reduce  the  effects  of  this  trade-off,  we  used  high  damping  both  for  anisotropy 
and  phase  velocities.  We  tested  our  inversion  results  by  using  a  variety  of  damping 
values.  However,  the  final  choice  of  damping  parameter  was  made  considering  model 
uncertainties  (Figure  20)  and  consistency  with  ANT  results.  Moreover,  the  features  we 
discuss  here  are  consistently  observed  in  all  inversions  regardless  of  damping. 
Additionally,  to  increase  reliability,  we  only  interpret  phase  velocities  with  errors  of 
about  0.06  km/s,  where  a  is  the  standard  deviation. 

The  third  and  last  step  of  the  inversion  is  obtaining  shear  wave  velocities  (Figure  16) 
using  the  anisotropic  phase  maps  from  the  second  step.  For  shear  wave  inversions,  we 
used  the  partial  derivatives  from  Saito  [1988].  For  the  initial  model,  the  AK135  model 
was  adapted  to  our  region  by  taking  average  Moho  depths  of  45  and  60  km  for  Qaidam 
Basin  and  the  rest  of  the  region  (Zhu  and  Helmberger,  1998),  respectively. 

We  defined  50  layers  for  shear  wave  inversion  from  the  surface  to  400  km  depth. 
According  to  Li  (2011),  who  used  a  similar  frequency  range  with  less  data,  the  shear 
wave  structure  cannot  be  determined  precisely  at  all  layers,  but  the  average  velocity  for 
each  defined  layer  can  be  resolved  for  depths  less  than  250  km.  Therefore,  we  do  not 
interpret  our  results  for  depths  >  200  km. 

Figures  15-16  show  our  tomographic  models  for  Rayleigh  wave  phase  and  shear  wave 
velocities.  Starting  from  the  south,  we  observe  a  high  velocity  body  (HVBl)  along  the 
BNS  and  further  south,  between  the  depths  of  ~60  and  170  km.  This  feature  is  highly 
heterogeneous  in  the  east-west  direction  with  a  shear  wave  anomaly  of  >  2%. 
Furthermore,  the  phase  velocity  maps  (Figure  15)  show  regions  of  relatively  low 
velocities  within  the  HVBl.  The  shear  wave  anomaly  cross-sections  show  that  the  HVBl 
has  a  sub-horizontal  geometry  dipping  northwards  slightly  beneath  the  BNS,  between 
~94E  and  97E.  A  further  isolated  high  velocity  body  (HVB2)  is  observed  at  ~91E-33N 
which  appears  to  be  connected  to  the  HVBl  beneath  the  BNS  at  ~ 130- 150  km,  becoming 
isolated  again  at  greater  depths.  This  vertically  continuous  feature  starts  at  ~60  km  depth 
(50  s)  and  gradually  increases  in  diameter  to  depths  greater  than  200  km  (Figures  15  and 
16). 

Further  north,  low  velocities  characterize  the  uppermost  mantle  beneath  the  N.  Qiangtang 
and  Songpan-Ganzi.  The  shear  wave  velocities  are  ~4%  slower  than  those  observed 
beneath  the  Lhasa  block  and  southern  Qiangtang  terrane  (Figure  14).  Anisotropic  fast 
directions  across  the  northern  and  southern  branches  of  the  Kunlun  Fault  tend  to  conform 
to  strikes  of  the  active  faults.  We  also  observe  the  largest  amplitude  of  azimuthal 
anisotropy  across  these  shear  zones.  Seismic  anisotropy  appears  to  be  consistent  at  depth 
with  the  exception  of  the  southeast  portion  of  our  region.  Similar  to  shear  wave  splitting 
measurements  (Huang  et  al.,  2000;  Le6n  Soto  et  al.,  2012),  azimuthal  fast  directions 
exhibit  a  general  clockwise  rotation  (Figure  15). 
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Period(s) 

Figure  17.  Standard  Error  (left  vertical  axis)  vs.  period  for  different  damping  parameter 
values.  Color  codes  are  shown  in  the  figure  legend.  Note  that  the  damping  and  damping 
parameter  are  inversely  proportional  (i.e.  lower  damping  parameter  values  mean  higher 
damping).  As  expected,  model  misfits  decrease  as  damping  increases  (damping 
parameter  decreases).  Considering  the  consistency  between  ANT  and  TPWT  results,  and 
the  model  misfits,  we  determined  that  0.2  is  the  most  suitable  damping  parameter  value. 
Moreover,  all  features  discussed  in  the  manuscript  are  consistently  observed  in  all 
inversions. 

Both  phase  maps  and  S  wave  anomalies  (Figures  15  and  16)  show  ~1%  higher  velocities, 
indicating  thicker  (~ 13 0-1 40  km)  and  possibly  rigid  lithosphere  beneath  the  Qaidam 
Basin.  We  also  observe  the  eastern  edge  of  the  Tarim  Basin.  In  approximate  agreement 
with  the  lithosphere-asthenosphere  boundary  (LAB)  phases  of  S-wave  receiver  functions 
[Zhao  et  al.,  2011],  the  Tarim  Basin  has  a  lithospheric  thickness  of  ~200  km,  while  the 
Qaidam  Basin  is  substantially  thinner.  Beneath  the  Kunlun  Shan,  south  of  Qaidam  Basin, 
we  observe  low  velocity  zones  that  are  concentrated  along  the  northern  and  southern 
branches  of  the  Kunlun  Fault,  extending  deeper  than  200  km. 

Comparisons  with  independent  phase  maps  using  ambient  noise  (Yang  et  al.,  2010) 
suggest  that  we  are  able  to  measure  the  phase  velocities  with  a  high  degree  of  precision 
(differences  less  than  0.06  km/s  within  the  array.  Figure  12).  Furthermore,  unlike 
previous  surface  and  body  wave  measurements  (Li  et  al.,  2008),  our  shear  wave 
observations  are  broadly  consistent  with  those  obtained  using  body  wave  tomography 
(Liang  et  al.,  2012).  Moreover,  the  comparison  of  the  splitting  data  from  Le6n  Soto  et  al. 
(2012)  and  our  azimuthal  fast  directions  for  143  s  shows  that  average  differences  are  ~1 
degree  and  10  degrees  for  Qiangtang  and  Songpan-Ganzi,  respectively  (Figure  15). 
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4.2  Pn  and  Sn  Travel  Time  Tomography 

We  created  a  reliable  tomographic  model  for  the  Pn  velocity  structure  in  western  China 
including  station  corrections  or  crustal  delay  times  (Figures  18  and  19).  Data  from  recent 
INDEPTH  IV  and  other  broadband  deployments  in  Tibet  are  used  to  invert  for  crustal 
delays  and  Pn  velocities  beneath  the  eastern  Tibetan  Plateau  and  surrounding  regions. 
Several  similar  high  velocity  regions  occur  north  of  the  eastern  Himalayan  syntaxis 
beneath  the  Lhasa,  Qiangtang,  and  Songpan-Ganzi  terranes.  We  interpret  these  as 
underthusted  Indian  continental  lithosphere  segments.  The  east-west  variability  of  Pn 
velocity  beneath  the  Himalayas  and  southern  Tibet  indicates  that  the  underthrusted  Indian 
continental  lithosphere  is  not  a  homogeneous  body.  Qaidam  Basin  along  with  the  Gonghe 
and  Xining  Basins  east  of  it  show  an  east-west  high  velocity  anomaly.  The  Indian 
lithosphere  extends  roughly  half-way  into  the  Tibetan  plateau,  with  more  penetration 
north  of  the  eastern  syntaxis. 
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Figure  18.  Our  final  tomographic  Pn  velocity  model  for  ray  paths  between  3  and  15 
degrees  using  data  from  the  Chinese  National  catalog  as  well  as  picks  from  all  temporary 
stations  deployed  throughout  the  Tibetan  plateau. 


Our  previous  Pn  tomography  studies  showed  a  contrast  between  high  velocities  beneath 
southern  Tibet  and  low  velocities  beneath  northern  Tibet,  but  we  now  find  that  this  model 
is  oversimplified.  Longer  ray  paths  show  a  vastly  different  velocity  structure  than  do  the 
shorter  ray  paths;  thus  there  is  significant  vertical  variation  in  the  velocity  structure  of  the 
mantle  lid.  Only  longer  ray  paths  show  low  velocity  beneath  northern  Songpan-Ganzi 
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terranes,  but  are  not  coincident  with  the  Bangong-Nujiang  suture  that  separates  northern 
and  southern  Tibet.  There  is  also  low  velocity  mantle  beneath  the  northeast  Songpan- 
Ganzi  Terrane  (in  Gansu  province).  Shorter  ray  paths  image  the  top  of  the  mantle  lid  and 
show  two  distinct  high  velocity  zones  beneath  the  Himalayas. 
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Figure  19.  Our  final  list  of  average  station  delays  for  all  stations  used  to  construct  the 
tomographic  model  shown  in  Figure  18. 
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Figure  20.  Our  final  tomographic  Pn  velocity  model  for  ray  paths  between  3  and  9 
degrees. 
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Figure  21.  Our  final  tomographic  Pn  velocity  model  for  ray  paths  between  6  and  12 
degrees. 
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Figure  22.  Our  final  tomographic  Pn  velocity  model  for  ray  paths  between  9  and  15 
degrees. 

In  addition  to  mapping  variations  in  Pn  velocity  we  have  also  created  a  model  for  Sn 
velocities  across  much  of  China  and  eastern  Tibet.  Furthermore,  in  order  to  be  able  to 
constrain  Vp/Vs,  a  key  parameter  for  earthquake  location,  for  this  region  we  have 
extracted  only  those  earthquake-station  pairs  with  both  Pn  and  Sn  travel  time  picks.  This 
will  allow  us  to  directly  compute  Vp/Vs  because  the  resolution  of  the  two  tomographic 
models  should  be  very  similar.  In  order  to  test  this  further,  however,  velocity  resolution 
tests  were  performed  to  determine  the  size  of  features  that  can  be  reliably  imaged,  and 
therefore  interpreted.  Checkerboard  tests  were  conducted  to  evaluate  the  effects  of  ray 
coverage  on  spatial  resolution.  A  test  checkerboard  velocity  model  was  created  by 
assigning  sinusoidal  velocity  anomalies  to  the  cells  of  the  model  domain  simultaneously. 
Synthetic  arrival  times  were  calculated  for  the  test  model  under  different  checkerboard 
sizes  with  the  same  number  of  earthquakes,  stations,  and  ray  paths  that  are  used  for  actual 
tomographic  inversion.  Gaussian  noise  was  added  to  the  synthetic  travel  times  and 
Gaussian  distributed  errors  with  a  standard  deviation  of  15  km  were  added  to  the  event 
locations  to  simulate  location  errors..  We  performed  tests  using  different  cell  sizes  for  the 
sinusoidal  velocity  anomalies.  The  spatial  resolution  is  considered  to  be  good  for  the 
regions  in  which  the  checkerboard  pattern  is  recovered.  The  tests  indicate  that,  for  most 
of  the  study  region,  velocity  anomalies  with  a  size  of  2.7°  x  2.7°  can  be  well  resolved, 
and  at  the  western  part  of  the  study  area,  the  spatial  resolution  can  reach  2.0°  x  2.0°  or 
even  better. 
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Figure  23.  A  map  of  eastern  Tibet  showing  the  major  physiographic  features  (left)  and  a 
plot  of  the  ray  path  coverage  where  we  observe  both  Pn  and  Sn  arrival  times  (Lii  et  al, 
2012). 


For  the  inversion  for  VpA^ s  for  the  uppermost  mantle,  we  used  the  Pn  and  Sn  phase  data 
for  epicentral  distances  between  3°  and  16°.  We  show  the  ray  path  coverage  in  Figure  23. 
The  initial  models  were  obtained  by  a  fit  to  the  travel  time-distance  curve  considering  the 
effect  of  the  global  curve.  The  average  Pn  and  Sn  velocities  of  the  upper  layer  are  8.10 
km/s  and  4.60  km/s  and  those  of  the  deeper  layer  are  8.19  km/s  and  4.65  km/s, 
respectively  (Figures  24  and  25).  The  travel  time  residuals  (difference  between  predicted 
and  reported  travel  times)  used  in  the  study  were  limited  to  a  range  of  ±5  s. 

In  general,  the  Pn  and  Sn  velocity  variation  patterns  are  similar  to  those  obtained  in  the 
model  of  the  upper  layer  discussed  above.  High  velocity  values  are  observed  beneath  the 
Sichuan  Basin  and  the  Indian  plate;  low  velocities  are  found  to  be  distributed  under  the 
Myanmar-Yunnan  region,  Hainan  region,  and  the  Songpan-Ganzi  fold  belt.  Four 
independent  datasets  of  the  inversion  of  the  Pn  and  Sn  travel  time  data  yielded  similar 
results  that  are  consistent  with  each  other,  and  this  confirm  the  credibility  of  the  inversion 
results..  In  our  model,  there  is  a  clear  velocity  difference  between  the  southern  and  the 
northern  parts  of  26°  N  in  the  Yunnan  region.  Low  velocities  are  found  south  of  26°  N. 
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Figure  24.  Sn  lateral  velocity  variations  for  the  uppermost  mantle  layer, 
average  Sn  velocity  is  4.60  km/s. 
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Figure  25.  Pn  lateral  velocity  variations  for  the  uppermost  mantle  layer.  Average 
Pn  velocity  is  8.10  km/s. 
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Figure  26.  A  tomographic  model  of  the  Vp/Vs  ratio  for  the  uppermost  mantle  beneath 
eastern  Tibet  and  the  surrounding  regions.  Average  velocity  ratio  is  1. 76.  Red  and  blue 
areas  correspond  to  ratios  higher  and  lower  than  the  average,  respectively. 


Figure  26  shows  the  VpA^s  ratio  distribution.  The  average  ratio  is  approximately  1.76 
beneath  the  study  area.  The  velocity  ratio  distributions  are  generally  similar  for  both  the 
Pn  and  the  Sn  models.  The  velocity  ratio  is  low  in  the  Sichuan  Basin.  In  the  Myanmar- 
Yunnan  and  Hainan  regions,  the  velocity  ratio  is  considerably  high,  indicating  that  the 
uppermost  mantle  had  a  high  temperature.  In  the  Yunnan  region,  there  are  some 
differences  in  the  velocity  ratio  distributions  between  the  two  layers.  From  the  heat  flow 
observation  (Wang  and  Huang,  1990;  Hu  et  ah,  2001)  and  cmst  velocity  structure  studies 
(Gao  et  ah,  2009;  Wu  and  Zhang,  2012;  Xu  et  ah,  2012),  we  infer  that  the  hot  material 
upwelling  is  not  vertical  but  is  more  complex  beneath  the  Yunnan  region.  For  the  upper 
layer,  Pn  and  Sn  were  determined  using  by  rays  with  epicenter  distances  between  3°  and 
10°.  The  average  Pn  and  Sn  velocities  of  this  layer  are  8.10  km/s  and  4.60  km/s, 
respectively.  The  average  VpWs  ratio  is  approximately  1.76  beneath  the  study  area.  The 
RMS  residuals  of  the  Pn  travel  times  decreased  from  2.0  s  to  1.1  s  after  inversion,  and  the 
residuals  of  Sn  travel  times  decreased  from  2.2  s  to  1.3  s.  In  general,  the  variations  of  the 
Pn  and  Sn  velocities  are  similar.  High  velocity  values  are  observed  beneath  the  Sichuan 


Approved  for  public  release;  distribution  is  unlimited. 


39 


Basin  and  the  Indian  plate.  Clear  low  velocities  are  found  under  the  Myanmar- 
Yunnan  region,  the  Hainan  region,  and  the  Songpan-Ganzi  fold  belt. 

At  the  Indian  plate  at  the  west  of  the  eastern  Himalayan  syntaxis  and  the  Sichuan  Basin, 
the  Pn  velocity  is  clearly  high  for  both  layers.  The  model  shows  that  the  lithospheres  of 
these  regions  are  cold  and  stable.  The  stable  lithospheres  of  the  Sichuan  Basin  block  the 
material  flow  from  the  Tibetan  Plateau  to  the  east  and  change  the  flow  direction  to 
the  south.  In  the  Myanmar- Yunnan  region,  low  Pn  velocities  are  distributed  around 
some  volcanoes  (Lei  et  al.,  2009a).  These  results  serve  as  further  evidence  of  the 
subduction  of  the  Indian  plate  beneath  the  Myanmar-Yunnan  region  (Li  et  al.,  2008), 
which  caused  the  back-arc  upwelling.  In  the  Hainan  region,  mantle  plumes  were  found 
using  body-wave  tomography  (Lei  et  al.,  2009b).  Our  inversion  models  also  support  this 
discovery.  The  Pn  velocity  under  the  Songpan-Ganzi  fold  belt  is  significantly  low, 
which  has  also  been  observed  in  previous  studies  (Liang  et  al.,  2004;  Pei  et  al.,  2007). 
This  low  velocity  can  be  primarily  attributed  to  high  temperatures,  and  secondarily  to 
compositional  variations  due  to  water  (Hearn  et  al.,  2004). 

In  general,  the  Pn  and  Sn  velocity  variation  patterns  are  similar  to  those  obtained  in  the 
model  of  the  upper  layer  discussed  above.  High  velocity  values  are  observed  beneath  the 
Sichuan  Basin  and  the  Indian  plate;  low  velocities  are  found  to  be  distributed  under  the 
Myanmar-Yunnan  region,  Hainan  region,  and  the  Songpan-Ganzi  fold  belt.  Four 
independent  datasets  of  the  inversion  of  the  Pn  and  Sn  travel  time  data  yielded  similar 
results  that  are  consistent  with  each  other,  and  this  confirms  the  credibility  of  the 
inversion  results.  A  previous  study  involving  an  SKS  analysis  revealed  that  the  fast 
direction  of  the  SKS  north  of  26°  N  is  N-S  and  south  of  26°  N  is  W-E  (Sol  et  al.,  2007). 
In  our  model,  there  is  a  clear  velocity  difference  south  and  north  of  26°  N  in  the  Yunnan 
region,  with  low  velocities  south  of  26°  N.  Considering  the  previous  SKS  and  Pn 
anisotropy  results  (Cui  and  Pei,  2009),  we  infer  that  the  regions  south  and  north  of  26°  N 
in  the  Yunnan  region  are  controlled  by  different  dynamic  processes.  The  southern  part  is 
affected  by  the  eastward  subduction  of  the  Indian  plate;  the  northern  part  is  affected  by 
material  flow  from  the  Tibetan  Plateau  (Zhang  et  al.,  2010;  Bai  et  al.,  2010). 
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Figure  27.  Map  showing  seismic  stations,  regional  seismic  events,  and  two-station  or 
two-event  paths  in  this  study.  Yellow  inverse  triangles  represent  the  Namche  Barwa 
experiment;  yellow  squares  represent  the  MIT-China  experiment;  red  triangles  represent 
the  INDEPTH-IV- ASCENT  experiment;  red  inverse  triangles  represent  the  INDEPTH- 
IV-UK  experiment;  red  stars  represent  the  NETS  experiment;  blue  diamonds  represent 
the  permanent  network  IC;  green  circles  represent  events.  The  cyan  and  black  lines 
represent  RTS  and  RTE  paths  at  1  Hz,  respectively.  The  epicenter  of  the  Wenchuan 
earthquake  is  marked  by  black  star,  and  the  Longmenshan  fault  (LMSF)  area  is  noted. 
Elevation  is  shown  with  a  gray  scale. 


4.3  Lg  and  Pg  O  Tomography 

We  used  more  than  16,000  Lg  spectra  combined  to  form  12,000  Reverse  Two  Station 
(RTM)  paths  that  cover  most  of  the  eastern  Tibetan  plateau.  We  analyzed  Lg  spectra 
between  0.5  Hz  and  2  Hz,  which  typically  represents  the  predominant  frequency  range 
for  Lg  phases  travelling  through  tectonically  active  crust.  We  have  selected  16  different 
frequency  bands  within  this  predominant  frequency  range  using  a  perturbation  of  0. 1  Hz, 
a  starting  central  frequency  of  0.5  Hz,  and  an  ending  central  frequency  of  2.0  Hz.  We 
used  a  bandwidth  of  0.8  Hz  for  each  of  these  frequency  bands.  A  moving  average  is  used 
to  smooth  each  spectrum,  and  a  cubic  spline  interpolation  is  used  to  replace  unstable  data 
points  if  at  least  80%  of  the  number  of  data  points  for  a  spectrum  fall  within  an  order  of 
magnitude  of  the  smoothed  spectrum. 
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The  study  area  is  divided  into  square  two  dimensional  eells  with  a  size  of  0.5°  x  0.5°  for 
our  tomographie  model.  A  set  of  maps  of  path  density  (i.e.  hit  counts)  have  been  built  for 
all  16  frequency  bands.  The  densest  path  coverage  is  found  within  an  area  where  the 
Namche  Barwa,  western  portion  of  ASCENT,  INDEPTH-IV-UK,  and  NETS  were 
deployed,  which  is  also  consistent  with  the  retrievable  area  of  2°  x  2°  anomalies  in  our 
resolution  tests  (Figure  10).  In  terms  of  the  tectonics,  this  area  mainly  consists  of  the 
Qaidam  Basin  (QB),  Songpan-Ganzi  terrane  (SGT,  or  Bayan  Har  block  [e.g.,  Vergne  et 
ah,  2002]),  Qiangtang  terrane  (QT),  Lhasa  terrane  (LT),  and  Eastern  Himalayan  Syntaxis 
(EHS)  (Figure  33). 
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Figure  28.  QLg  (a)  and  Qpg  (b)  tomographic  images  at  1  Hz.  Thin  solid  lines  represent 
faults,  and  thick  lines  represent  tectonic  sutures.  ANHF:  Anninghe  Fault;  BNS:  Bangong- 
Nujiang  Suture;  EHS:  Eastern  Himalayan  Syntaxis;  IP:  Indian  Plate;  lYS:  Indus-Yalu 
Suture;  JS:  Jinsha  Suture;  KL:  Kunlun  mountain  belt;  LB:  Lhasa  Block;  LMS: 
Longmenshan  mountain  belt;  MJF:  Minjiang  Fault;  QB:  Qaidam  Basin;  QLGB:  Qinghai 
Lake-Gonghe  Basin;  QL:  Qinling  mountain  belt;  QLS-NS:  Qilian  Shan-Nan  Shan 
Mountain  belt;  QT:  Qiangtang  Terrane;  SB:  Sichuan  Basin;  SGFB:  Songpan-Ganzi  fold 
belt;  XSHF:  Xianshuihe  Fault;  YGF:  Yushu-Ganzi  Fault. 
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Figure  29.  Models  of  QLg  at  0.5  Hz  (a),  1  Hz  (b),  and  2  Hz  (c). 


A  high  QLgo  zone  is  loeated  near  the  Eastern  Himalayan  Syntaxis  (EHS)  and  probably 
extends  eastward  to  the  Bangong-Nujiang  Suture  (BNS)  in  the  southeastern  TP  (Figure 
28),  whieh  was  also  observed  by  Li  et  al.,  2008.  The  Kunlun  mountain  range  is  at  the 
northern  margin  of  the  TP,  where  we  observe  a  nearly  eontinuous  E-W  very  low  QLgo 
“band”  that  divides  the  high  QLgo  Qaidam  Basin  and  low  to  middle  QLg  central  TP. 
Another  low  QLg  zone,  nearly  along  the  western  Longmenshan  thrust  belt,  laterally 
divides  the  high  QLg  Sichuan  Basin  and  the  low  to  middle  QLg  eastern  TP. 
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Theoretical  studies  suggest  that  a  complex  velocity  structure  leads  to  scattering 
attenuation  and  regional  phase  mode  conversion  when  Lg  propagates  in  a  strongly 
heterogeneous  crust  (e.g.,  Kennett,  1989).  But  practically,  large  numbers  of  studies  of  Lg 
propagation  imply  that  the  QLg  is  a  good  approach  for  crustal  shear  wave  Qs,  and  that 
scattering  plays  only  a  minor  role  and  is  probably  negligible  at  low  frequencies 
(Herrmann  and  Kijko,  1983;  Mitchell,  1995;  Sarker  and  Abers,  1998;  Jemberie  and 
Nyblade,  2009).  We  only  see  Lg  frequencies  lower  than  10  Hz  in  our  study,  which 
suggests  that  the  conditions  of  measuring  QLg  and  its  frequency  dependence  are  similar  to 
the  studies  referenced  above.  Meanwhile,  focusing  and  defocusing,  especially  lateral 
refraction,  probably  affects  the  fundamental  assumptions  of  the  RTM  and  the  accuracy  of 
tomographic  modeling  (Mitchell,  1995).  However,  the  mechanism  of  focusing  and 
defocusing  is  mostly  caused  by  small-scale  velocity  anomalies,  probably  less  than  10  km 
(Sarker  and  Abers,  1998),  which  is  secondary  and  probably  negligible  for  measurements 
from  relatively  longer  paths  such  as  our  minimum  inter-station  distance  of  150  km.  Thus 
the  lateral  variation  of  Qlso  probably  represents  the  lateral  variation  of  crustal  Qs,  and 
furthermore  contains  information  about  geothermal  structures  and  active  tectonics  in  the 
crust  (e.g.  Mitchell,  1995). 

Our  observations  of  high  QLgo  zones,  such  as  the  Qaidam  Basin,  probably  represent 
tectonically  relatively  stable  crust  with  low  temperatures.  The  very  low  QLgo  band  along 
the  Kunlun  mountain  range  dramatically  corresponds  with  Cenozoic  volcanic  fields 
developed  during  Indo-Asia  Collision  (Yin,  2010)  and  regions  with  high  horizontal  shear 
strain  rate  modeled  from  Quaternary  fault  slip  and  GPS  measurements  (Holt  et  ah,  2000); 
this  correlation  suggests  a  relationship  among  high  shear  strain  rate,  young  volcanism, 
and  Q  in  the  crust.  Unlike  prior  studies  noting  that  the  Qiangtang  crust  perhaps  contains 
partial  melting  in  northern  Tibet,  our  model  suggests  that  the  high  shear  strain 
accumulated  in  the  crust  of  the  Songpan-Ganzi  fold  belt  and  Kunlun  leads  to  high 
temperature,  and  possibly  partially  molten  crust. 
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Figure  30.  Lg  site  term  maps  at  frequencies  of  0.5  Hz  (a),  1  Hz  (b),  and  2  Hz  (c).  We 
observe  significant  site  term  differences  between  the  Tibetan  Plateau  and  major  basins, 
probably  dependent  on  topography  or  sedimentary  thickness. 


We  also  observe  that  the  site  responses  have  a  very  strong  spatial  structure  across  the  TP; 
we  consistently  observe  high  site  amplification  values  within  the  Qaidam  Basin  and  low 
amplification  values  within  high  TP  (Figure  30).  This  does  not  correlate  with  any 
difference  in  instrument  response  or  suggest  a  potential  artificial  effect  from  deployment 
style.  There  is  no  contradiction  between  the  pattern  of  lateral  variation  of  measured  site 
response  in  this  study  and  the  assumption  of  the  frequency  dependency  of  site  responses 
(e.g.,  Huang  et  al.,  2005)  because  we  calculate  relative  site  terms.  But  such  a  regular 
pattern  is  unlikely  consistent  with  the  assumption  that  the  site  responses  are  strongly 
dominated  by  near-surface  shear  wave  velocity  and  seismic  structures  (Rodgers  et  al., 
2006),  and  its  lateral  variation  does  not  correspond  with  the  distribution  of  different 
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geological  strata.  In  contrast,  this  pattern  shows  strong  evidence  that  the  site 
amplification  is  dependent  on  topography  (Zhang  and  Lay,  1994)  or  sedimentary 
thickness  (Mitchell  and  Hwang,  1987).  This  suggestion  also  helps  to  explain  the  different 
QLgo  model  measured  from  TSM,  in  which  the  site  amplification  is  neglected  and  the  QLgo 
measurements  are  biased.  In  these  models,  QLgo  within  the  Qaidam  Basin  is  not  as 
distinguishable  from  the  Lg  attenuation  in  the  high  TP. 

In  terms  of  our  source  terms,  we  observe  a  relatively  linear  relationship  between 
estimated  logarithmic  RTM  source  terms  and  earthquake  magnitude  Ml  (Figure  31)  when 
3<Mi<7.  This  shows  a  possible  application  of  RTM  Lg  amplitude  in  accurately  detecting 
earthquake  magnitude,  because  in  theory  it  eliminates  the  contribution  from  geometrical 
spreading,  site  response,  focusing/defocusing,  scattering,  and  attenuation  in  propagation 
within  inhomogeneous  crust.  We  emphasize  that  its  practical  application  is  probably  only 
limited  by  seismic  station  distribution. 

Similar  to  QLg  as  a  proxy  of  crustal  Qs,  the  Qpg  of  another  crustal  regional  phase  Pg  is 
assumed  as  a  proxy  of  crustal  Qp.  Although  the  propagation  of  Pg  still  contains 
uncertainties  in  theory  (Kennett,  1989),  many  studies  have  suggested  that  it  is  reasonable 
and  satisfactory  to  consider  Pg  phases  as  post-critical  P  waves  trapped  in  a  crustal  wave 
guide  (e.g.  Hearn  and  Clayton,  1986;  Steck  et  al.,  2009).  We  apply  a  very  similar  data 
processing  flow  to  the  Pg  data  as  we  do  with  the  Lg  data;  while  some  parameters,  such  as 
the  group  velocity  window,  the  geometrical  spreading  coefficient,  the  threshold  of 
maximum  inter-station  distance,  etc,  have  been  revised  to  appropriate  values  for  Pg  (Bao 
et  al.,  2011).  The  calculated  RTM  Qpg  tomographic  image  is  strongly  correlated  with  the 
RTM  QLg  model  (Figure  28b):  high  attenuation  regions  throughout  most  of  the  northern 
TP  and  low  attenuation  regions  surrounding  TP.  However,  we  observe  differences  in  the 
relative  attenuation  structure  between  Pg  and  Lg  in  the  Qilian  Shan-Nan  Shan  ranges  and 
central  Qiangtang  Terrane  that  probably  indicate  rheological  anomalies  in  the  crust. 
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Figure  31.  The  linear  relationship  between  normalized  logarithmic  source  terms  solved 
from  Lg  RTM  and  magnitude  Ml.  The  solid  line  represents  the  least-square  linear  fit  line 
and  the  two  dashed  lines  parallel  to  it  represent  its  95%  confidence  interval.  It  indicates 
that  the  source  term  solved  from  Lg  amplitude  using  RTM  is  well  correlated  with  seismic 
magnitude  Ml  when  3  ^Mi  ^7.  This  method  eliminates  the  contribution  of  path,  site,  and 
source  radiation  pattern,  and  its  practical  application  is  only  limited  by  seismic  station 
alignment. 
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Figure  32  Examples  of  the  RTM  HQig-cp  correlation.  The  x-axis  represents  path 
azimuth  cp  and  y-axis  represents  l/Qig  value.  The  values  of  l/Qig  are  normalized  using 
integration  within  an  interval  of  3.6°  (1/50  of  180°).  In  each  solvable  interval,  the 
normalized  l/Qig  is  assumed  as  the  median  number  (black  dot),  and  error  bar  is  also 
shown  as  gray.  The  black  curves  represent  fitting  curves  using  sinusoidal  function  of  2(p 
with  a  corresponding  phase,  whose  period  is  180°.  The  latitude  and  longitude  of  each 
block  are  specified  in  the  title. 
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Figure  33.  Imaged  variations  of  isotropic  Q  {\IA),  anisotropy  magnitude  (5),  and  high-g 
direction  in  the  crust  of  eastern  Tibetan  Plateau  at  different  frequencies  of  (a)  0.6  Hz, 
(b)  0.8  Hz,  (c)  1  Hz,  and  (d)  1.2  Hz.  The  color  scale  represents  the  variations  of  1/A,  and 
the  length  and  orientation  of  black  vectors  represent  1/B  and  6,  respectively.  Major  faults 
and  sutures  are  also  shown  using  the  database  HimaTibetMap-1 .0  (Taylor  and  Yin,  2009), 
while  strike-slip  faults  are  green,  other  types  of  faults  are  gray,  and  sutures  are  yellow. 
QB:  Qaidam  basin;  SGFB:  Songpan-Ganzi  fold  belt;  QT:  Qiangtang  terrane;  LB:  Lhasa 
block;  EHS:  Eastern  Himalyan  Syntaxis. 
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Figure  34.  Comparison  of  azimuthal  anisotropy  of  l/Qig,  for  different  frequencies,  in 
the  eastern  Tibetan  Plateau.  The  fast  directions  are  overlain  at  each  grid  point,  with 
vector  length  scaled  by  magnitude  of  the  anisotropy,  for  16  narrow  frequency  bands 
between  0.5  Hz  and  2.0  Hz  as  described  in  the  text.  At  most  points  all  the  vectors  are 
essentially  parallel,  indicating  very  weak  frequency  dependency  within  most  of  this 
region. 

Our  images  showing  lateral  variations  in  A,  B,  and  0  at  discrete  frequencies  from  0.5  to 
2.0  Hz,  in  discrete  steps  of  0.1  Hz,  are  solved  from  the  l/Qig  azimuthal  anisotropy 
tomography.  Fig.  36  shows  four  examples  of  these  tomographic  images  at  different 
frequencies  of  0.6  Hz,  0.8  Hz,  1  Hz,  and  1.2  Hz.  Only  the  cells  with  more  than  5  hit- 
counts  are  included.  To  compare  the  results  with  tectonics,  we  employ  the 
HimaTibetMap-1.0  database  [Taylor  and  Yin,  2009]  to  show  major  faults  and  plate 
boundaries.  Three  aspects  of  the  model,  on  the  anisotropic  part  (containing  two  variables 
of  B  and  0)  and  isotropic  part  (A)  of  l/Qig,  are  discussed  and  interpreted. 
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(i)  The  major  pattern  of  anisotropie  magnitude  (5)  is  that  there  is  a  strong  azimuthal 
anisotropy  of  MQig  within  the  northern  TP.  The  maximum  B  of  the  whole  study  region, 
typically  higher  than  0.01,  is  observed  within  the  northwestern  SGT,  the  Kunlun  Shan 
(mountain  range),  and  northern  QT.  Meanwhile,  B  is  mostly  less  than  0.004  in  the  central 
QB,  which  shows  a  very  significant  discontinuity  between  the  northern  and  southern  side 
of  the  Kunlun  Shan.  B  slowly  decreases  to  below  0.005  in  the  central  and  southern  TP, 
and  appears  to  increase  slightly  around  the  EHS.  In  general  we  observe  a  very  weak 
frequency-independent  B  (Fig.  8).  (ii)  A  nearly  uniformly  W-E  high-Q  direction  (^  is 
observed  in  the  northern  TP.  The  highest  B  values  are  found  in  the  northwestern  SGT, 
with  ^correlating  strongly  with  the  strike  of  the  major  transform  fault  planes,  such  as  the 
North  Kunlun  Fault  (NKF)  and  South  Kunlun  Fault  (SKF).  The  central  and  southern  TP 
are  probably  divided  into  three  sections  in  terms  of  the  distribution  of  6.  A  clockwise 
rotation  of  0  is  around  the  EHS,  with  a  northwest-southeast  (NW-SE)  direction  on  the 
western  side  of  EHS  and  a  nearly  N-S  direction  on  the  eastern  side  of  EHS.  We  also 
observe  that  0  is  largely  frequency  independent  (Fig.  8).  (iii)  A  strong  frequency- 
dependent  A  in  the  eastern  TP  is  solved  by  isolating  the  anisotropic  part  of  HQig,  which 
leads  to  the  major  differences  among  Fig.  36  (a),  (b),  (c),  and  (d).  At  0.6  Hz,  the  central 
QB  has  HA  higher  than  400  and  the  northwestern  SGT,  northern  QT,  and  central  LT  have 
MA  probably  lower  than  150.  The  1/^4  values  at  0.8  Hz  are  typically  higher  than  those  at 
0.6  Hz  with  a  shift  of  50.  For  frequencies  larger  than  1  Hz,  we  find  three  relatively  high 
HA  zones  (350  at  1  Hz  and  >400  at  1.2  Hz)  in  the  southern  QT,  EHS,  and  the  center  of 
northwestern  SGT. 

There  are  likely  two  possible  interpretations  for  the  strong  azimuthal  dependence  of  HQig 
along  the  Kunlun  Shan:  the  influence  from  non-uniform  geometry  of  the  crust,  or  from 
the  crustal  anisotropy.  Surface  topography  [Wu  and  Wu,  2001]  and  Moho  geometry, 
sediment  thickness  [e.g.,  Kennett,  1986],  and  inhomogeneous  velocity  structures  of  the 
crust  can  generate  a  complicated  seismic  wave  field  at  high  frequencies  (>  0.5  Hz). 
Across  the  Kunlun  Shan,  steep  topography  (>  2  km),  a  thick  sedimentary  basin  (>  5  km 
in  the  Qaidam),  and  rapid  Moho  depth  changes  (15-20  km)  [e.g.,  Zhu  and  Helmberger, 
1998;  Karplus  et  al,  2011]  all  could  lead  to  substantial  scattering  of  the  Eg  wave  field. 
Previous  studies  have  observed  an  energy  loss  of  Lg  when  its  path  laterally  crosses  the 
Kunlun  Shan  [Zhao  et  al.,  2003],  which  is  probably  associated  with  large  geometrical 
fluctuations  along  this  mountain  range.  A  relationship  between  ray  path  azimuth  a  and 
the  energy  loss  of  Lg  has  been  suggested  by  focused  and  defocused  synthetic  ray  paths 
[Bostock  and  Kennett,  1990].  But  the  role  of  the  influence  of  geometry  on  seismic 
attenuation  from  a  generalized  focusing  and  defocusing  is  actually  still  undeterminable 
for  Lg  [Furumura  and  Kennett,  1997].  The  size  of  attenuation  anomalies  in  our 
tomographic  model  is  approximately  200  km  and  the  focusing  and  defocusing  produced 
by  substantial  seismic  velocity  variations  should  probably  create  smaller  scale  Q 
anomalies  [Barker  and  Abers,  1998].  The  site  responses  beneath  all  seismic  stations 
solved  by  the  RTM,  not  included  in  the  attenuation  term,  has  been  suggested  to  represent 
the  influences  from  topography  or  sediment  thickness  [Bao  et  a/.,  2011].  Our  model  with 
high  anisotropy  in  other  parts  of  the  TP  does  not  clearly  show  any  strong  correlation  with 
Moho  geometry  [e.g.,  Wittlinger  et  al,  1996].  We  suggest  that  generalized  focusing  and 
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defocusing  does  not  provide  a  necessary  and  sufficient  condition  for  the  azimuthal 
anisotropy  of  HQig- 

The  seismic  attenuation  measured  for  the  upper  crust  has  been  also  suggested  to  include 
the  effect  from  faults  and  cracks  associated  with  fluids  {Mitchell,  1995].  Fault  zones  have 
been  suggested  to  play  an  important  role  in  attenuation  of  coda  waves  resulting  from 
scattering  {Jin  et  al,  1994],  Although  it  is  still  unknown  how  Lg  samples  the  crust  at 
different  depths,  the  Qig  value  should  be  affected  by  a  crustal-scale  fault  zone.  An 
azimuthal  anisotropy  of  llQig  is  expected  to  be  observed  if  such  a  fault  system  contains 
well  aligned  fractures  and  interstitial  fluids.  There  is  no  evidence  of  abundant  fracture 
fluids  beneath  TP,  where  a  very  dry  crust  is  probably  widespread  {Hacker  et  al,  2000]; 
however,  such  an  interpretation  with  fracture  fluids  is  unnecessary.  It  has  been  presented 
that  the  azimuthal  anisotropy  of  HQig  can  be  described  as  transversely  isotropic  media 
(TIM)  containing  multiple  layers  or  aligned  fractures  with  horizontal  symmetry  axis  {Zhu 
and  Tsvanskin,  2006;  Chichinina  et  al,  2009],  and  therefore  be  caused  by  the  layers  and 
faults  themselves,  not  the  content  of  fluid  within  them.  Its  anisotropic  magnitude  is 
related  to,  but  much  stronger  than,  the  magnitude  of  the  velocity  anisotropy  {Liu  et  al, 
2005;  Chichinina  et  al,  2006;  Chapman,  2009].  This  type  of  attenuation  from  multiple 
layers  and  aligned  fractures  is  actually  a  shape  preferred  orientation,  and  is  one  of  the 
major  mechanisms  for  the  crustal  anisotropy  {Werner  and  Shapiro,  1998]. 

The  NKF  and  SKF  are  the  two  eastward  extended  branches  of  the  Kunlun  Fault  (KF)  and 
both  are  major  left-lateral  strike-slip  faults  in  the  northeastern  TP.  After  separating  at 
approximately  94°E,  the  NKF  nearly  follows  the  E-W  orientation  of  the  KF  while  the 
SKF  turns  to  southeast  and  nearly  parallels  the  orientation  of  the  Jinsha  Suture  (JS),  the 
tectonic  boundary  of  the  QT  and  SGT.  We  observe  an  alignment  of  the  High  Q  directions 
along  both  the  southern  and  northern  branches  of  the  Kunlun  fault  zones:  E-W  along  the 
northern  branch  and  NW-SE  along  the  southern  branch.  The  tectonic  evolution  of  the 
SGT  crust  has  likely  been  very  complicated  since  the  Mesozoic  after  the  closing  of  the 
Songpan  Ocean,  including  major  amounts  of  crustal  shortening  and  thickening  {Roger  et 
al,  2010].  The  complicated  distribution  of  fractures  and  faults  within  the  folded  crust  of 
the  SGT  may  play  an  important  role  in  the  generation  of  the  observed  anisotropic  velocity 
and  attenuation  structures.  Both  the  NKF  and  SKF  likely  extend  into  the  heterogeneous 
upper  crust  probably  with  very  low  velocity  and  Q  between  them  {Wang  et  al,  2009; 
Karplus  et  al,  2011].  If  the  crustal  anisotropy  is  the  major  source  of  the  azimuthal 
anisotropy  of  MQig  estimated  in  our  model,  the  upper  crust  of  the  northwestern  comer  of 
the  SGT  may  be  assumed  nearly  transversely  isotropic  with  abundant  parallel  faults  and 
fractures  following  the  orientations  of  the  NKF,  SKF,  and  JS.  Such  a  pattern  is  probably 
generated  by  high  rates  of  deformation  within  the  SGT  that  have  been  measured  by 
observations  of  Quaternary  fault  slip  rates  and  GPS  observations  [e.g..  Holt  et  al,  2000]. 
We  consistently  observe  strong  azimuthal  anisotropy  of  CQig  with  consistently  east- west 
directions  in  the  northern  portion  of  the  QT.  This  orientation  is  consistent  with  the 
characteristics  of  upper  cmstal  anisotropy  within  the  SGT  {Vergne  et  al,  2003]  and  about 
3°  west  of  our  model  {Ozacar  and  Zandt,  2004],  which  suggests  a  common  origin  for 
both  the  velocity  and  attenuation. 
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The  pattern  of  high-g  direetions  0  is  dramatieally  similar  with  the  fast  wave  direction  in 
the  upper  to  middle  crustal  models  (about  20-30  km  in  our  discussion)  of  azimuthal 
anisotropy  of  Rayleigh  waves  [5^  et  ah,  2008;  Yi  et  al,  2010;  Yao  et  al,  2010].  These 
models  have  different  lateral  resolutions  (5°x5°  in  Su  et  al,  2008,  4°x4°  in  Yi  et  al,  2010, 
and  2°x2°  in  Yao  et  al,  2010)  not  higher  than  our  model  (2°x2°).  They  show 
corresponding  magnitudes  of  velocity  anisotropy  of  approximately  3%-5%  and  fast  wave 
directions  of  nearly  W-E  in  the  western  portion  with  a  clockwise  rotation  in  the  eastern 
portion  in  our  study  area  at  an  upper  to  middle  crustal  depth.  A  more  complicated  pattern 
of  fast  wave  directions  around  the  EHS,  like  the  pattern  of  our  high-g  directions,  is 
observed  from  the  model  of  Yao  et  al,  2010  with  the  highest  lateral  resolution  in  the  three 
models  of  Rayleigh  waves.  The  similarity  between  the  models  of  azimuthal  anisotropy  of 
HQig  and  shear  wave  velocity  further  suggests  a  probable  pattern  of  crustal  anisotropy 
and  supports  the  assumed  relationship  between  velocity  anisotropy  and  attenuation 
anisotropy.  If  the  estimated  high-g  direction  correlates  with  the  orientation  of  crustal- 
scale  fault  systems,  its  lateral  pattern,  including  W-E  direction  in  the  northern  TP,  the  N-S 
direction  to  the  north  of  the  EHS,  and  the  clockwise  rotation  around  the  EHS,  it  may 
represent  the  pattern  of  upper  and  likely  middle  crustal  fabric  caused  by  faults  and 
fractures,  which  is  probably  related  to  shape-preferred  orientation  and  different  major 
deformation  motions  in  the  whole  TP.  The  nearly  W-E  trend  of  high-Q  direction  0  in  the 
northeastern  TP  is  also  close  to  the  patterns  of  Global  Position  System  velocity  vector 
measurements  [e.g.  Wang  et  al,  2001]  and  shear  wave  splitting  fast  directions  [e.g., 
Huang  et  al,  2000;  Wang  et  al,  2008;  Li  et  al,  2011].  This  correlation  supports  the 
assumption  of  mechanical  coupling  between  the  crust  and  lithospheric  mantle  during 
their  deformation  in  the  northeastern  TP.  But  the  southern  TP  shows  a  complexity  where 
the  stability  of  such  coupling  is  likely  affected  by  the  probably  underthrusting  Indian 
lithosphere.  Meanwhile,  since  the  propagation  of  Lg  waves  is  also  likely  affected  by  the 
complicated  lower  crustal  velocity  structures  within  the  southern  TP  [Ndbelek  et  al, 
2009],  a  possible  contribution  from  the  lattice  preferred  orientation  may  lead  to  a  more 
complicated  situation  in  the  representation  of  the  azimuthal  anisotropy  of  MQig  for  the 
estimated  crustal  anisotropy. 

4.4  Sn  Q  Tomography  (Tibet  and  Northeastern  China) 

We  studied  Sn  propagation  across  the  Tibetan  Plateau  and  northeastern  China  in  order  to 
quantify  upper  mantle  shear  wave  attenuation.  Efficient  and  inefficient  longer  period  Sn 
waves  were  identified  using  data  from  the  INDEPTH-IV  array.  More  than  12,000 
waveforms  were  reviewed  and  the  Sn  propagation  efficiencies  determined  from 
instruments  used  in  the  ASCENT  and  NETS  arrays  in  the  northeastern  Tibetan  plateau 
and  a  further  6,000  waveforms  were  reviewed  for  our  study  in  northeastern  China.  We 
were  able  to  identify  more  than  2000  efficient  paths  from  data  that  were  bandpass  filtered 
through  0.2  and  0.05  Hz  for  long  period  Sn  waveforms.  We  used  both  the  two  and 
reverse  two  station  method  to  calculate  Sn  Q  that  eliminates  any  source  contribution  to 
the  Sn  waveform.  All  waveforms  were  also  corrected  for  instrument  response.  This 
method  allowed  us  to  tomographically  determine  variations  in  Sn  Q  across  the  Tibetan 
plateau. 
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Figure  35.  Sn  efficiency  map  for  Tibet  with  all  efficient  paths  shown  on  top  of  the 
blocked  and  inefficient  paths. 

Our  results  suggest  that  Sn  Q  is  low  in  and  around  the  active  fault  systems  like  the 
Kunlun  Fault.  This  observation  is  also  consistent  with  shear-wave  velocity  models  for  the 
uppermost  mantle.  This  observation  points  to  a  thin  and  hot  lithosphere  in  northeastern 
Tibet  up  to  the  Qaidam  basin.  Both  tomographic  RTM  and  TSM  models  show  portions  of 
high  Q  (low  attenuation)  zones  in  central  Tibet  which  are  interpreted  to  be  the  leading 
edge  of  the  underthrusting  Indian  lithopshere.  Within  the  Qaidam  basin  we  observe  two 
regions  of  very  high  Sn  Q  that  suggest  that  we  are  either  smearing  the  effects  from  the 
Qilian  Shan  hot  uppermost  mantle  or  there  are  two  distinct  domains  of  the  Qaidam 
mantle. 
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Figure  36.  Sn  Q  tomography  of  northeastern  Tibet  at  1  Hz.  We  once  again 
observe  the  uppermost  mantle  beneath  the  Songpan  Ganzi  to  be  highly 
attenuation  which  is  largely  consistent  with  P-wave  and  S-wave  velocity 
estimates  of  the  same  region. 

Comparisons  of  the  RTM  and  TSM  (Figures  36  and  37)  models  for  Sn  Q  across  eastern 
Tibet  show  significant  differences.  This  observation  would  suggest  that  the  crustal  legs 
of  Sn  are  having  a  significant  impact  on  measured  Sn  amplitudes. 

In  addition  to  working  on  long  period  Sn  Q  studies  in  Tibet,  we  extended  our  work  to 
include  northern  China  with  the  eventual  goal  in  future  work  of  creating  an  Sn  Q  map  for 
all  of  China  and  some  of  the  surrounding  regions.  All  of  the  seismograms  in  our  study 
have  been  manually  reviewed  by  applying  two  bandpass  filters  with  comer  frequencies  of 
0.5  to  2  Hz  and  0.1  to  0.5  Hz  to  identify  Sn.  We  classified  each  seismogram  into  one  of 
three  categories:  no  Sn  (blocked),  inefficient  Sn  and  efficient  Sn.  If  the  Sn  wavetrain  can 
be  seen  for  both  high  and  frequencies  seismograms,  we  categorize  it  as  efficient.  If  no 
distinguishable  high-frequency  Sn  wavetrain  is  observed  but  Sn  could  be  observed  in  the 
low  frequency  band,  or  an  ambiguous  signal  was  found,  we  classified  it  as  inefficient. 
Otherwise,  we  assigned  the  seismogram  as  “no  or  blocked  Sn”.  Figure  38  shows  all  Sn 
ray  paths  with  color  coded  efficiencies  used  in  our  study.  It  shows  that  the  inefficient  Sn 
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ray  paths  are  primarily  restricted  to  those  paths  that  cross  the  Sea  of  Japan,  and  the 
efficient  Sn  ray  paths  are  mostly  found  for  continental  paths  in  our  study  region. 
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Figure  37.  Top:  RTM  Sn  Q  tomography  of  northeastern  Tibet  at  1  Hz;  Bottom; 
RTM  Sn  Q  tomography  of  northeastern  Tibet  at  0.5  Hz.  Note  the  differences 
between  the  Sn  TSM  and  RTM  methods,  which  suggest  that  the  station  and  event 
site  amplification  have  a  significant  impact  on  Sn  amplitudes. 
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Figure  38:  A  map  showing  the  Sn  raypaths  across  northeastern  China  using  broadband 
data  from  the  NECESS  array.  The  propagation  efficiencies  are  labeled  by  color:  Blue, 
green  and  red  lines  present  efficient,  inefficient  and  “no  or  blocked  Sn  ”  raypaths, 
respectively. 

We  inspected  individual  seismograms  to  remove  the  data  that  are  too  noisy  or  have 
incorrect  timing.  The  Pn  arrival  time  was  automatically  calculated  and  we  also  manually 
re-picked  Pn  for  a  more  accurate  arrival  time.  It  was  used  to  obtain  a  pre-Pn  window  and 
calculate  the  ambient  noise  level.  Only  those  estimated  Sn  spectra  with  signal-to-noise 
ratio  larger  than  2  were  kept.  Furthermore,  we  cut  a  pre-Sn  window  between  4.8  and  5.0 
km/s  and  an  extensional  Sn  window  between  4.1  and  4.8  km/s.  If  the  root  mean  square  of 
the  extensional  Sn  signal  to  that  of  the  pre-Sn  was  less  than  1 .05,  the  data  were  removed. 
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Figure  39.  Sample  record  section  for  an  earthquake  occuring  along  the  Pacific 
subduction  zone. 


Approved  for  public  release;  distribution  is  unlimited. 


59 


The  Sn  velocity  window  was  defined  for  4.3  to  4.7  km/s  uppermost  mantle  velocities. 
The  Sn  window  was  automatically  estimated  to  start  with  a  fast  velocity  of  4.8  km/s  and 
end  with  a  slow  velocity  of  4. 1  km/s.  To  better  isolate  the  Sn  spectra,  we  adjusted  the  Sn 
windows  manually.  A  window  with  a  Pg  velocity  of  5.1  km/s  and  Lg  group  velocity  of 
3.5  km/s  was  used  to  avoid  mistaking  Pg  coda  or  early  Lg  for  Sn.  Figure  7  shows  some  of 
the  efficient  Sn  seismograms  for  one  event.  We  can  see  that  the  Sn  arrival  time  is  around 
the  center  of  the  Sn  window,  which  indicates  the  Sn  velocity  is  approximately  4.45  km/s. 
In  our  study,  we  set  an  Sn  velocity  of  4.4  km/s  based  on  the  velocity  model  provided  by 
Zheng  et  al.  (2011). 


45‘ 


45’ 


The  calculated  Qo  values  were  used  to  map  the  lateral  variation  of  Qo  by  applying  the 
LSQR  algorithm  (Figure  42).  We  can  roughly  divide  the  Qo  values  into  four  regions:  low 
(Qo  <400),  medium  low  (400<  Qo<800),  medium  high  (800<  Qo<1200),  and  high 
(Qo>1400).  We  observed  low  Qo  in  the  Xilinhote  block,  the  Songnen  block,  the  eastern 
Zhangguangcailing  block,  and  the  Yanji  block  which  are  part  of  the  Songnen  unit, 
medium  low  Qo  in  the  Erguna-Xing’an  block,  medium  high  Qo  in  the  central  and  northern 
Songliao  basin,  and  high  Qo  in  the  southern  Songliao  basin,  Jiamusi  block  and  part  of 
Harlar-Tamsag  basin  and  Erlian  basin.  In  general  we  observe  that  the  basins  are 
characterized  by  high  Sn  Qo  values.  The  borders  between  high  Qo  and  low  Qo  are  mostly 
along  the  active  faults. 
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Figure  41  An  example  of  linear  regression  fitting  for  interstation  RTM  Qo  values  for  the  same 
station  pair  NE79-NEA2.  The  title  of  each  plot  indicates  the  station  pair,  event  pair,  QO  values 
and  errors.  Note  the  large  differences  in  the  slopes  indicating  the  relatively  poor  control  we  have 
on  the  frequency  dependence  for  this  particular  two  station  pair.  We  can,  however,  use  these 
types  of  repeated  paths  to  estimate  confidence  limits  on  our  Q  models. 
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We  measured  the  RTM  Sn  Q  values  at  discrete  frequencies  of  0.5  Hz,  1  Hz  and  2  Hz.  The 
data  procedure  is  revised  from  Xie  (2002a)  and  Xie  et  al.  (2004). ).  Furthermore  we  did 
not  obtain  any  effective  RTEM  paths  due  to  the  distribution  of  the  stations  and  events. 

We  set  the  minimum  inter-station  distance  of  150  km  in  order  to  minimize  Q  estimation 
errors.  Also,  we  set  the  maximum  azimuthal  angle  at  15°  to  minimize  modeling  errors 
and  the  effects  of  multi-pathing.  We  set  the  minimum  and  maximum  allowable 
interstation  Sn  Q  values  to  10  and  3000  respectively.  Only  Q  values  with  error  less  than 
50%  and  that  met  the  interstation  Qo  quality  control  criteria  were  accepted.  In  the  end  we 
obtained  8003  reverse  two-station  paths  (Figure  40)  and  Qo  values.  The  Qo  and  q  values 
and  deviations  were  calculated  from  a  linear  regression.  Examples  of  interstation  spectral 
ratios  and  the  Qo  values  are  shown  in  figure  41.  We  found  that  the  Qo  values  range  from 
277  to  576. 
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Figure  42.  RTM  Sn  Q  map  at  1  Hz  with  a  cell  size  of0.5°x0.5°.  Blue  triangles  are 
stations. 


The  lateral  variation  of  Sn  Qo  was  helpful  to  infer  the  origin  of  seismic  anomalies  in  the 
upper  mantle.  By  correlating  our  Sn  Q  model  with  velocity  measurements,  we  could 
determine  which  originated  from  temperature  and  which  originated  from  compositional 
variations  in  the  uppermost  mantle.  Liu  et  al.  (2001)  illustrated  that  the  pattern  of 


Approved  for  public  release;  distribution  is  unlimited. 


62 


magmatism  since  4.5  Ma  would  suggest  that  the  mantle  upwelling  strongly  affected 
attenuation  of  Sn  in  regions  within  the  Songliao  basin. 

We  also  measured  TSM  Q  in  this  region  in  order  to  eompare  it  to  our  RTM  Q  model. 
Any  differenees  would  likely  originate  from  the  site  terms  listed  in  equation  (19).  When 
preparing  the  two-station  paths,  we  kept  those  paths  with  the  more  distant  station  having 
an  Sn  effieiency  of  “no  Sn”  and  obtained  1685  two-station  paths  (Figure  40).  We 
calculated  TSM  Sn  Qo  for  0. 5-1.5  Hz  (Figure  41).  Compared  with  our  RTM  model, 
generally,  it  shows  higher  Qo  values  in  high  Q  regions  and  lower  Qo  values  in  low  Q 
regions.  This  may  be  due  to  the  site  effeet  that  was  not  eliminated  in  the  TSM  Sn  Qo. 
Typieally,  high  site  responses  occur  in  basins,  which  may  lead  to  higher  TSM  Qo  values 
in  the  basin. 


Figure  43.  Cheekerboard  test  with  a  2°x2°  cells  using  3000  reverse  two-station  paths. 
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Figure  44.  Sn  Q  at  (a)  0.5  Hz  and  (b)  2.0  Hz  with  O.S'^xQ.S'^  cells.  Blue  triangles 
denote  stations. 
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A  2D  checkerboard  test  with  10%  Gaussian  distributed  noise  shows  relatively  good 
resolution  for  a  2°  x  2°  anomaly  size  in  most  of  the  study  area  (Figure  43).  In  general 
we  do  observe  some  east  west  smearing  because  of  the  predominance  of  east-west 
reverse  two  station  paths.  Despite  this,  we  are  still  able  to  resolve  Q  anomalies  at  a  2 
degree  cell  size. 

We  have  broadened  the  frequency  band  to  0.1~2.0  Hz  for  measuring  q.  We  found  low  q 
values  in  the  center  of  our  study  area  where  we  suggest  that  anelastic  attenuation 
dominates,  whereas  high  rj  values  around  the  edge  suggest  scattering  attenuation  or  more 
likely  worse  resolution  for  the  frequency  dependence.  At  the  edges  of  our  model  we  have 
fewer  crossing  paths  and  which  leads  to  larger  errors  and  possibly  larger  values  of  q. 
Stronger  frequency  dependence  would  be  more  consistent  with  scattering  attenuation 
rather  than  intrinsic  attenuation.  In  general  we  observe  only  a  weak  frequency 
dependence  which  surprisingly  suggests  that  intrinsic  attenuation  dominates  throughout 
northeastern  China. 


The  frequency  dependence  of  Sn  Q  appears  to  be  quite  complicated.  In  general  we  do 
observe  a  strong  correlation  between  the  frequency  dependence  and  the  Q  value  itself, 
which  is  commonly  observed  for  high  frequency  regional  phase  Q.  This  observation 
suggests  that  at  least  part  of  the  frequency  dependence  is  a  systematic  problem  with 
measuring  larger  Q  values.  We  can  see  this  effect  in  the  Songliao  Basin  where  we  see  the 
same  spatial  patterns  in  q  as  we  do  in  Q 
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Figure  45.  Our  percentage  error  of  Sn  Q  at  1  Hz  across  our  study  region.  It  is  necessary 
to  show  the  Q  error  as  a  percentage  because  larger  Q  values  are  significantly  less  well 
constrained  than  lower  Q  values. 
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We  have  found  a  very  large  number  of  reverse  two  station  paths  (i.e.  more  than  3000), 
however,  many  of  these  paths  are  repeated  paths.  In  faet  we  have  approximately  800 
unique  reverse  two  station  paths  for  the  NECESS  array  (mostly  east-west).  This  means 
we  have  a  large  number  of  repeat  paths  for  eaeh  of  our  reverse-two  station  pairs.  An 
example  of  these  repeated  paths  is  shown  in  Figure  41  for  station  pair  NE79-NEA2.  We 
ean  use  these  types  of  repeated  paths  to  estimate  confidenee  limits  of  our  regional  phase 
Q  models  by  measuring  the  standard  deviations  of  these  repeated  paths.  We  have 
eompiled  standard  deviations  for  all  two  station  pairs  with  at  least  10  measurements  and 
then  eomputed  a  standard  deviation  for  each  pair.  We  then  tomographically  mapped  the 
spatial  variation  in  errors  to  produce  error  maps.  Figures  45  and  46  show  our  model  error 
estimates  for  the  Q  at  1  Hz  (Qo)  and  the  frequency  dependence  (q)  respectively.  The 
errors  for  the  frequency  dependence  show  particularly  large  values  in  regions  with  high  Q 
which  again  suggests  that  we  have  a  difficult  time  constraining  the  frequency  dependence 
for  large  Q  values,  as  was  seen  in  figure  41.  In  regions  with  lower  Q  values,  on  either 
side  of  the  Songliao  Basin,  we  do  appear  to  be  able  to  constrain  the  frequency 
dependence  reasonably  well. 
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Figure  46.  The  error  for  the  Q  frequency  dependence  term  q  across  our  study  region. 
Note  the  very  large  error  along  the  Songliao  Basin  which  suggests  that  we  have  little  to 
no  control  on  Q’s  frequency  dependence  at  Q  values  larger  than  1000 for  Sn  amplitudes. 
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4.5  ML  Amplitude  Tomography 


The  inversion  eontained  station  gains,  event  gain  eorrections,  geometric  spreading,  and 
regionally  varying  Q.  While  both  the  amplitude  and  measured  period  were  included  in 
the  inversion,  the  period  data  often  were  reported  with  only  one  significant  digit  of 
precision.  Since  ML  periods  are  small  (0. 1  to  2  seconds)  and  in  the  denominator,  this 
means  that  errors  in  the  period  could  cause  large  errors  in  the  predicted  amplitude.  We 
attempted  to  use  it  in  the  inversion  but  results  were  unreliable.  As  a  result,  we  fixed  the 
periods  to  0.5  seconds  in  the  inversion.  Nevertheless,  it  is  important  to  note  that 
assumptions  about  frequency  dependence  greatly  affect  both  the  geometric  spreading  and 
attenuation  values. 


Distance  {krr) 


Figure  47.  Log  amplitude  versus  offset  for  the  raw  data  set  before  selection.  Amplitudes 
were  first  converted  to  units  of  nanometers  and  then  normalized  for  event  size  by 
subtracting  the  magnitude.  Our  fit  to  these  data  found  a  spreading  coefficient  of  0. 75 
and  an  average  Q  of 490.  This  absolute  Q  value  could  be  biased  by  the  assumed  velocity 
and  period;  however,  the  pattern  of  attenuation  found  in  the  tomography  is  stable. 


The  inversion  restricted  data  to  events  with  depths  less  than  33  km  and  source-receiver 
distances  between  0.5  and  9  degrees.  Data  had  an  average  period  of  0.5  seconds  so  we 
restricted  the  data  to  those  between  0. 1  and  1 .0  seconds.  In  addition,  we  iteratively 
eliminated  stations  with  less  than  10  arrivals  and  events  with  less  than  25  recording 
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stations  and  also  iteratively  eliminated  arrivals  with  residuals  greater  than  one  order  of 
magnitude.  We  have  assumed  an  average  wave  speed  of  3  km/s  and  the  average  period 
of  0.5  seconds. 


The  results  gave  a  geometrical  spreading  factor  of  0.8  and  an  average  Q  value  of  484 
(Figure  47).  The  estimated  Q  value  depends  inversely  on  the  assumed  wave  velocity  and 
period.  We  attempted  to  use  the  measured  period  in  the  inversion,  but  the  results  were 
unreasonable  with  high  spreading  and  Q  values.  This  occurs  because  the  period  data 
often  were  reported  with  only  one  significant  digit  of  precision  and  are  subject  to  errors. 
Furthermore,  since  ML  periods  are  small  (0.1  to  21  seconds)  and  in  the  denominator, 
small  errors  in  the  period  can  result  in  large  errors  in  the  estimated  amplitudes.  Attempts 
at  including  a  frequency  dependent  Q  were  also  made;  however  these  resulted  in  the 
untenable  frequency  dependences.  Again,  we  believe  this  is  because  of  the  poor 
precision  and  accuracy  of  the  period  measurement.  We  did  find,  however,  that  including 
frequency  dependent  Q  could  significantly  alter  the  spreading  and  attenuation  estimates. 
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Figure  48.  Tomographic  image  of  Q  across  Tibet  using  amplitudes  from  the  Chinese 
Earthquake  Bulletin. 


A  comparison  with  our  work  using  Lg  and  Pg  data  is  somewhat  limited  because  of  the 
somewhat  limited  overlap  between  the  two  models,  however,  we  do  observe  consistency 
across  the  Qaidam  basin  (high  Q)  and  Qilian  Shan  mountain  belt  (low  Q).  Furthermore, 
we  have  found  similar  patterns  to  the  relative  station  terms  between  the  Ml  tomography 
and  RTM  approaches.  We  consistently  observe  a  de-amplification  of  crustal  phases  in 
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the  Tibetan  plateau  while  we  observe  amplification  in  the  Qaidam  basin  and  surrounding 
regions  (Figure  49).  The  absolute  amplitude  of  the  Q  anomalies  are  somewhat  different 
however.  The  Ml  amplitude  tomography  (shown  in  Figure  48)  had  larger  Q  values 
(Q~1000)  than  we  found  using  the  RTM  and  TSM  approaches  (Q~500). 
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Figure  49.  Relative  station  terms  for  the  tomography.  Station  terms  are  the  logarithm  of 
the  station  gains  and  are  expressed  in  magnitude  units. 

4.6  Joint  Inversion  of  Receiver  Functions  and  Surface  Wave  Dispersion  Curves 
We  conducted  a  feasibility  study  of  using  our  3-D  shear  wave  velocity  model  for  Tibet  by 
doing  a  joint  inversion  of  dispersion  curve  data  using  both  phase  and  group  velocity  data 
combined  with  published  receiver  functions  across  the  ASCENT  seismic  network.  We 
compared  different  velocity  models  using  independent  techniques  to  see  which  portions 
of  the  models  are  common  to  all  of  our  models.  We  used  the  receiver  function  work  of 
Han  et  al.  (2012)  and  surface  wave  phase  velocities  of  Ceylan  et  al.  (2013)  combined 
with  the  global  group  velocities  from  Pasyanos  (2010)  by  applying  the  joint  inversion 
method  of  Julia,  et  al.  (2000).  We  used  Rayleigh  wave  phase  velocities  obtained  from  the 
two  plane  wave  method  and  receiver  functions  calculated  using  the  method  of  Liggoria 
and  Ammon,  1999.  We  took  advantage  of  the  constraints  provided  by  PS  move-out  in 
receiver  functions  by  creating  a  large  number  of  receiver  function  stacks  over  a  relatively 
narrow  distance  range  (~5  degrees)  and  then  simultaneously  modeled  all  of  the  stacks. 
By  utilizing  the  PS  move-out  this  way  we  could  better  constrain  the  Vp/Vs  ratio  for  the 
crust.  We  were  generally  able  to  fit  both  data  reasonably  well  with  the  exception  of  some 
of  the  PS  converted  phase  amplitudes.  An  example  of  the  data  fits  is  shown  in  Figure  50 
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for  a  station  located  directly  on  the  Kunlun  fault  zone  just  south  of  Golmud  China.  This 
result  is  typical  of  what  we  have  found  for  these  stations  in  this  region. 

We  have  obtained  72  one  dimensional  velocity  models  for  all  of  the  ASCENT  array 
stations;  interpolating  between  these  72  models  we  were  able  to  construct  a  three 
dimensional  velocity  model.  Furthermore,  we  estimated  a  crustal  thickness  value  for  all 
72  stations  by  choosing  the  center  of  the  large  positive  velocity  gradient  at  the  base  of  the 
crust.  We  found  that  by  choosing  the  center  of  the  velocity  gradient  we  obtained  results 
roughly  consistent  with  Han  et  ah,  2012.  If  we  chose  the  base  or  beginning  we  obtained 
generally  unrealistic  crustal  thickness  values.  Figure  5 1  is  our  optimal  Moho  map  using 
the  center  depth  of  the  1-D  velocity  gradients;  this  model  shows  cmstal  thickness  patterns 
consistent  in  a  relative  sense  with  the  migrated  receiver  functions:  a  thinning  of  the  crust 
toward  the  Qaidam  basin  and  the  thickest  crust  south  of  the  Jinsha  suture.  We  also  found 
further  evidence  of  an  upper  crustal  low  velocity  zone  and  a  crustal  thickness  of 
approximately  70  km  within  the  Kunlun  Fault  zone.  We  observed  a  reasonable 
agreement  between  mantle  P-to-S  converted  waves  and  the  LVZ  imaged  by  the  surface 
wave  data. 
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station  COOl 


Figure  50.  The  results  of  one  of  our  Joint  Inversions  for  station  COOl.  This  station  was 
located  near  the  northern  edge  of  the  Tibetan  plateau  with  the  Qiadam  basin  located 
just  to  the  north.  We  found  evidence  of  a  low  velocity  layer  in  the  uppermost  mantle  that 
is  consistent  with  the  surface  tomography  of  Ceylan  et  al,  2012  and  our  Pn 
tomographic  models  using  shorter  ray  paths. 
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Figure  51.  Our  Moho  map  generated  from  the  optimal  1-D  velocity  models  from  70 
stations  across  northeastern  Tibet.  The  triangles  show  the  location  of  all  stations  were 
obtained  optima  1-D  models  for.  This  map  was  constructed  by  taking  the  middle  of  the 
steepest  velocity  gradient  from  the  lower  crust  to  the  uppermost  mantle. 
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5.  CONCLUSIONS 

5.1  Seismic  Velocity  Studies 


Figure  52.  A  rough  quantitative  comparison  of  the  overlapping  regions  of  the 
body  and  surface  wave  models  of  the  shear  wave  velocity.  We  normalized  the 
velocity  perturbations  by  the  maximum  percent  change  (negative  or  positive)  and 
then  subtracted  the  two  normalized  models  to  produce  the  contour  maps. 
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We  compared  different  velocity  models  using  different  techniques  to  see  which  portions 
of  the  models  are  robust.  We  created  three  dimensional  P  and  S  wave  velocity  structure 
models  for  eastern  Tibet  using  both  finite  frequency  body  wave  tomography  (not  on  this 
contract)  and  surface  wave  tomography  (funded  by  this  contract).  In  general  the  upper 
mantle,  from  the  Main  Frontal  Thrust  to  central  Tibet  around  33°-34°N,  exhibits  laterally 
variable  P-  and  S-wave  velocity  anomalies  extending  to  at  least  250  km  in  depth. 
Significant  alternating  low  and  high  velocity  anomalies  are  observed  beneath  southern 
Tibet,  elongated  in  a  north-south  direction  and  extending  to  at  least  150  km  depth, 
possibly  deeper.  These  studies  confirm  that  low-velocity  Tibetan  crust  ends  abruptly  at 
approximately  105  degrees  longitude.  These  low  velocity  structures  are  consistent  with 
measurements  of  low  Lg  and  Pg  Q  values,  indicating  that  these  anomalies  are  a  result  of 
high  temperatures  and  not  just  composition.  The  low  velocities  show  a  strong  correlation 
with  the  strain  rates  computed  from  GPS  and  Quaternary  fault  data  suggesting  that  strain 
heating  plays  some  role  in  the  generation  of  the  lowest  velocities.  It  is  also  important  to 
note  that  neither  the  attenuation  nor  the  velocity  structure  is  symmetric  with  respect  to  the 
major  sutures  within  Tibet,  and  overall  we  observe  a  much  more  complex  velocity 
structure  than  previously  thought. 

Overall  we  observe  a  fairly  strong  correlation  between  the  different  features  in  our 
various  seismic  models  for  the  eastern  half  of  the  Tibetan  plateau.  In  order  to  test  which 
parts  of  our  models  differ  significantly  we  have  attempted  to  do  a  direct  comparison 
between  the  finite  frequency  S-wave  travel  time  models  and  our  surface  wave  velocity 
models.  Direct  comparison  of  these  models  is  impossible  because  the  finite  frequency 
models  only  offer  relative  velocities  while  the  surface  wave  models  give  us  absolute 
velocity.  Therefore  we  converted  our  surface  wave  models  to  a  relative  model  and  then 
normalized  both  of  the  models  to  either  +1.0  or  -1.0%  depending  upon  which  is  larger. 
We  then  subtracted  the  two  normalized  relative  S-wave  velocity  models  and  the  results 
are  shown  in  Figure  52.  The  largest  differences  between  the  two  models  occur  at  the 
Bangong  Nujiang  suture  (BNS).  It  is  important  to  note  that  we  observe  a  steeply  dipping 
high  velocity  anomaly  in  both  body  and  surface  wave  model  but  the  location  of  the 
anomalies  is  different. 

Large  two  dimensional  arrays  in  eastern  Tibet  have  yielded  a  number  of  important  results 
and  the  proliferation  of  seismic  instrument  pools  in  China  offer  a  tremendous  opportunity 
to  further  improve  the  resolution  of  three  dimensional  models  within  China.  In  particular 
large  aperture,  close  spaced  arrays  of  broadband  stations  across  the  margins  of  the 
Tibetan  plateau  should  help  to  address  fundamental  questions  about  the  different  modes 
of  deformation  within  the  plateau. 

5.2  Interpretation  of  the  Pn  and  Sn  Tomography 

In  general,  the  Pn  and  Sn  velocity  variation  patterns  are  similar  to  those  obtained  in  the 
model  of  the  upper  layer  discussed  above.  High  velocity  values  are  observed  beneath  the 
Sichuan  Basin  and  the  Indian  plate;  low  velocities  are  distributed  under  the  Myanmar- 
Yunnan  region,  Hainan  region,  and  the  Songpan-Ganzi  fold  belt.  Four  independent  data 
sets  of  the  inversion  of  the  Pn  and  Sn  travel  time  data  yielded  similar  results  that  are 
consistent  with  each  other,  and  this  confirms  the  credibility  of  the  inversion  results.  A 
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previous  study  involving  an  SKS  analysis  revealed  that  the  fast  direction  of  the  SKS 
north  of  26°  N  is  N-S  and  south  of  26°  N  is  W-E  (Sol  et  al.,  2007).  In  our  model,  there  is 
a  clear  velocity  difference  between  the  southern  and  the  northern  parts  of  26°  N  in  the 
Yunnan  region.  Low  velocities  are  found  south  of  26°  N.  Considering  the  previous  SKS 
and  Pn  anisotropy  results  (Cui  and  Pei,  2009),  we  infer  that  the  regions  north  and  south  of 
26°  N  in  the  Yunnan  region  are  controlled  by  different  dynamic  processes.  The  southern 
part  is  affected  by  the  eastward  ]  subduction  of  the  Indian  plate  and  the  northern  part  is 
affected  by  the  material  flow  from  the  Tibetan  Plateau  (Zhang  et  al.,  2011;  Bai  et  al., 
2010). 

In  eastern  Tibet  and  throughout  eastern  China  the  average  Vp/Vs  ratio  is  approximately 
1.76  beneath  the  study  area  (lower  than  the  global  average).  The  velocity  ratio 
distributions  are  generally  similar  for  the  Pn  and  the  Sn  models.  The  velocity  ratio  is  low 
in  the  Sichuan  Basin.  In  the  Myanmar- Yunnan  and  Hainan  regions,  the  velocity  ratio  is 
quite  high,  indicating  that  the  uppermost  mantle  has  a  high  temperature.  In  the  Yunnan 
region,  there  are  some  differences  in  the  velocity  ratio  distributions  between  the  two 
layers.  From  the  heat  flow  observations  (Wang  and  Huang,  1990;  Hu  et  al.,  2001)  and 
crustal  velocity  structure  studies  (Gao  et  al.,  2009;  Wu  and  Zhang,  2012;  Xu  et  al.,  2012), 
we  infer  that  the  hot  material  upwelling  is  not  vertical  but  is  more  complex  beneath  the 
Yunnan  region,  and  the  deeper  layer  is  significantly  affected  by  the  hot  material 
upwelling. 

In  the  back-arc  region  around  Myanmar,  significant  low  Pn  and  Sn  velocities  were  found 
to  be  distributed  around  volcanoes.  This  is  strong  evidence  that  the  lithosphere  of  the 
Indian  plate  subducted  into  the  mantle  beneath  Myanmar  and  caused  the  hot  material 
upwelling  at  the  back  arc.  The  high  temperature  or  even  partial  melting  in  the  uppermost 
mantle  caused  the  low  velocity  and  the  high  Vp/Vg  ratio.  Surprisingly  this  is  less  than  the 
global  average  for  the  uppermost  mantle  which  is  approximately  1.79.  This  number  is 
similar  to  Vp/Vs  ratios  that  are  found  in  many  subduction  zones,  suggesting  that  perhaps 
the  mineralogy  of  the  eastern  Tibetan  upper  mantle  is  a  result  of  the  subduction  of  the 
Tethyan  ocean  in  the  Mesozoic. 

The  uppermost  mantle  velocities  of  the  Sichuan  Basin  and  the  Songpan-Ganzi  fold  belt 
are  quite  different.  The  high  Pn  and  Sn  velocities  of  the  Sichuan  Basin  imply  that  the 
lithosphere  in  this  region  is  stable  and  cold,  and  thus  the  weaker  and  likely  more  ductile 
material  beneath  the  Songpan  Ganzi  terranes  flow  moving  to  the  south  along  the  edge  of 
the  strong  Sichuan  basin  lithospheric  mantle.  The  low  velocity  in  the  Songpan-Ganzi 
region  is  consistent  with  that  obtained  in  the  previous  studies.  The  very  low  Vp/Vs  ratios 
beneath  the  Sichuan  block  is  typical  for  many  cratonic  regions  and  can  be  attributed  in 
part  to  a  depleted  lithospheric  mantle. 


There  are  significant  velocity  differences  north  and  south  of  26°  N  in  the  Yunnan  region. 
From  the  SKS  anisotropy  result,  we  inferred  that  these  parts  were  controlled  by  different 
dynamic  processes.  The  southern  part  was  affected  by  the  eastward  subduction  of  the 
Indian  plate,  and  the  northern  part  was  affected  by  the  material  flow  from  the  Tibetan 
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Plateau.  The  Pn  and  Sn  velocities  are  low  and  the  Vp/Vs  ratio  is  high  in  the  Hainan 
region,  which  offers  further  evidence  of  mantle  plumes  beneath  the  Hainan  region.  The 
differences  in  the  velocity  ratio  distributions  between  the  models  of  the  two  layers  in  the 
Yunnan  region  indicates  that  the  hot  materiel  upwelling  was  not  vertical  but  is  more 
complex  and  the  deeper  layer  was  significantly  affected  by  the  hot  material  upwelling. 

5.3  0  Tomography 

Among  various  methods  used  to  measure  crustal  Q,  RTM  appears  to  be  the  most  reliable. 
Its  RTS  case  does  have  the  shortcoming  that  it  is  geometrically  limited  by  the  distribution 
of  seismic  stations,  but  this  is  improved  by  the  RTE  case,  which  makes  RTM  more 
practical  to  effectively  map  variations  crustal  Q  using  regional  phases.  The  estimation  of 
site  terms  using  RTM  in  this  study  indicates  that  the  site  response  is  probably  strongly 
dependent  on  topography  and/or  sedimentary  thickness.  The  linear  relationship  between 
estimated  logarithmic  source  terms  and  earthquake  magnitudes  Ml  when  shows 

that  RTM  Lg  amplitude  is  applicable  in  accurately  estimating  earthquake  magnitude.  We 
create  high-resolution  (for  1°  x  1°  to  2°  x  2°  anomalies)  QLg  and  Qpg  models  of  the 
eastern  TP  and  adjacent  areas  using  RTM  in  Q  measurements  and  the  LSQR  algorithm  in 
tomography.  The  QLg  and  Qpg  tomographic  models  show  strong  lateral  variation  in  the 
eastern  TP  and  adjacent  areas.  And  low  Q  zones  correlate  with  Cenozoic  volcanic  fields 
and  active  tectonics,  which  is  circumstantial  evidence  that  we  are  effectively  isolating  the 
path  based  attenuation  of  high  frequency  Lg  and  Pg  waves. 

Furthermore  we  measure  significant  azimuthal  anisotropy  of  1/QLg  within  the  eastern 
Tibetan  Plateau  using  a  reverse  two-station/event  method.  Attenuation  anisotropy  appears 
to  have  similar  orientations  to  the  azimuthal  velocity  anisotropy  measured  using  surface 
waves.  A  series  of  path-based  l/Qpg  values  are  used  to  frame  an  inverse  problem  in  order 
to  solve  for  the  frequency-dependent  isotropic  and  anisotropic  portions  of  l/Qpg- 
Resolution  tests  show  that  this  model  may  retrieve  anomalies  as  small  as  2°  x  2°. 
Anisotropic  l/Qpg  is  strong  in  the  northwestern  SGFB  with  the  high-Q  direction 
paralleling  major  strike-slip  fault  planes.  This  suggests  that  the  crustal  velocity  and 
attenuation  anisotropy  have  common  origins.  The  very  large  azimuthally  dependent 
fluctuations  in  l/Qpg  (Figure  35)  suggest  that  for  Lg  propagation  in  Tibet,  it  is  probably 
important  to  include  an  anisotropic  component  in  all  Q  models  in  order  to  better  predict 
high  frequency  amplitudes.  These  fluctuations  are  large  enough  to  block  Lg  for  some 
azimuths  and  propagate  inefficiently  for  different  path  directions  across  the  same  part  of 
Tibet.  We  have  also  assumed  that  the  isotropic  1/QLg  is  a  result  of  intrinsic  attenuation 
and  scattering  caused  by  randomly  distributed  scatterers  and  applied  a  model  of 
Dainty  (1981)  to  estimate  the  intrinsic  attenuation  of  QLg  which  is  highly  dependent 
on  the  temperature  of  the  crust.  The  crustal  temperature  estimated  for  the  Qiangtang 
terrane  in  this  study  corresponds  with  the  results  of  previous  studies.  The  northwestern 
comer  of  SGFB  may  have  a  high  cmstal  temperature  possibly  resulting  from  strain 
heating  of  this  strongly  deforming  flysch  complex. 

We  also  demonstrated  that  it  is  feasible  to  measure  high  frequency  Sn  Q  using  the  RTM 
method.  We  showed  that  at  least  for  Q  the  measurements  are  repeatable,  especially  for 
lower  values  of  Sn  Q  (see  Figure  48).  However,  for  measuring  the  frequency  dependence 
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our  results  are  not  as  robust,  suggesting  we  cannot  constrain  the  frequency  dependence 
for  Q  values  larger  than  1000.  Based  on  the  broadband  seismograms  from  regional  Sn 
waves,  we  determined  an  optimal  tomographic  Sn  Q  model  of  the  uppermost  mantle 
beneath  NE  China.  When  determining  the  Q  values  from  regional  phases,  it  is  essential  to 
eliminate  the  source  response,  site  effects,  radiation  pattern,  and  instrument  response. 
Although  the  TSM  does  not  effectively  remove  relative  site  responses  in  Q  values,  it  is 
widely  utilized  as  the  two  stations  give  it  a  decided  advantage  over  the  single  station 
methods  and  comparing  with  the  RTM,  it  has  fewer  assumptions  that  may  be  violated. 
The  RTM  computes  the  most  accurate  Q  values  since  it  removes  the  site  effect;  however, 
the  rigorous  requirements  for  applying  the  RTM  limits  the  scope  of  its  application  except 
when  we  have  seismic  networks  such  as  the  NECESS  array. 

The  frequency  dependence  of  Q  is  complicated  and  the  power-law  may  be  over 
simplified.  Within  a  narrow  frequency  band,  it  may  be  an  acceptable  approximation.  For 
our  study,  the  Sn  Q  appears  consistent  with  the  rj  tomography.  Negative  r]  are  shown  near 
the  high  7;  this  may  be  a  result  of  the  change  of  the  waveguide  shape  and  thus  Sn 
attenuation  is  dominated  by  scattering.  Our  Sn  Q  tomography  results  show  significant 
lateral  variation  in  NE  China,  which  may  reflect  the  rheology  in  the  uppermost  mantle. 
The  Sn  Q  values  are  relatively  high  in  the  Songliao  basin,  which  suggests  lower 
temperature  and  lower  viscosity  than  in  the  surrounding  area.  The  lower  Sn  Q  values 
correspond  with  Quaternary  volcanism.  Our  Sn  Q  model  is  consistent  with  previous  shear 
wave  velocity  models  at  the  Moho.  The  high  Q  regions  correspond  to  the  high  velocity 
area  and  the  low  Q  region  corresponds  to  the  low  velocity  area. 
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RTM 

Reverse  Two  Station  Method 

TSM 

Two  Station  Method 

TPWP 

Two  Plane  Wave 

TPWT 

Two  Plate  Wave  Tomography 

TP 

Tibetan  Plateau 

NKF 

Norhtem  Kunlun  Fault  zone 

SKF 

Southern  Kunlun  Fault  zone 

SGT 

Songpan  Ganzi  Terrane 

GPS 

Global  Positioning  System 

SKF 

Southern  Kunlun  Fault  zone 

QT 

Qiangtang  Terrane 

BNS 

Bangong  Nujiang  Suture 

JS 

Jinsha  Suture 

Approved  for  public  release;  distribution  is  unlimited. 


87 


DISTRIBUTION  LIST 


DTIC/OCP 

8725  John  J.  Kingman  Rd,  Suite  0944 
Ft  Belvoir,  VA  22060-6218  1  cy 

AFRL/RVIL 

Kirtland  AFB,  NM  87117-5776  2  cys 

Official  Record  Copy 

AFRL/RVBYE/Robert  Raistrick  1  cy 
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